diff --git a/src/TurbulenceConvection_deprecated/EDMF_functions.jl b/src/TurbulenceConvection_deprecated/EDMF_functions.jl index 9692e843cb1..41c40ca5c3a 100644 --- a/src/TurbulenceConvection_deprecated/EDMF_functions.jl +++ b/src/TurbulenceConvection_deprecated/EDMF_functions.jl @@ -97,20 +97,20 @@ function compute_sgs_flux!( @inbounds for i in 1:N_up a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) Ifau = CCO.InterpolateC2F(; a_up_bcs...) - a_up = aux_up[i].area - w_up_i = prog_up_f[i].w - @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up_i - w_gm) + a_up = aux_up.:($i).area + w_up_i = prog_up_f.:($i).w + @. aux_up_f.:($$i).massflux = ρ_f * Ifau(a_up) * (w_up_i - w_gm) @. massflux_h += ρ_f * ( Ifau(a_up) * (w_up_i - w_gm) * - (If(aux_up[i].h_tot) - If(h_tot_gm)) + (If(aux_up.:($$i).h_tot) - If(h_tot_gm)) ) @. massflux_qt += ρ_f * ( Ifau(a_up) * (w_up_i - w_gm) * - (If(aux_up[i].q_tot) - If(q_tot_gm)) + (If(aux_up.:($$i).q_tot) - If(q_tot_gm)) ) end @@ -253,7 +253,7 @@ function set_edmf_surface_bc( q_surf = q_surface_bc(grid, state, edmf, i, param_set) a_surf = area_surface_bc(surface_conditions, edmf) e_kin = @. LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), + C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f.:($$i).w))), ) / 2 e_pot_surf = geopotential(thermo_params, grid.zc.z[kc_surf]) ts_up_i_surf = @@ -264,10 +264,11 @@ function set_edmf_surface_bc( e_kin[kc_surf], e_pot_surf, ) - prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf - prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf - prog_up_f[i].w[kf_surf] = CCG.Covariant3Vector( + prog_up.:($i).ρarea[kc_surf] = ρ_c[kc_surf] * a_surf + prog_up.:($i).ρae_tot[kc_surf] = + prog_up.:($i).ρarea[kc_surf] * e_tot_surf + prog_up.:($i).ρaq_tot[kc_surf] = prog_up.:($i).ρarea[kc_surf] * q_surf + prog_up_f.:($i).w[kf_surf] = CCG.Covariant3Vector( CCG.WVector(FT(0)), CC.Fields.local_geometry_field(axes(prog_up_f))[kf_surf], ) @@ -388,21 +389,21 @@ function compute_implicit_up_tendencies!( LBF = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(CCG.WVector(FT(0)))) @inbounds for i in 1:N_up - w_up = prog_up_f[i].w + w_up = prog_up_f.:($i).w - ρarea = prog_up[i].ρarea - ρae_tot = prog_up[i].ρae_tot - ρaq_tot = prog_up[i].ρaq_tot + ρarea = prog_up.:($i).ρarea + ρae_tot = prog_up.:($i).ρae_tot + ρaq_tot = prog_up.:($i).ρaq_tot - tends_ρarea = tendencies_up[i].ρarea - tends_ρae_tot = tendencies_up[i].ρae_tot - tends_ρaq_tot = tendencies_up[i].ρaq_tot + tends_ρarea = tendencies_up.:($i).ρarea + tends_ρae_tot = tendencies_up.:($i).ρae_tot + tends_ρaq_tot = tendencies_up.:($i).ρaq_tot volume_term = @. -p_c / ρ_c * (-(∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea)))) @. tends_ρarea += -∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea)) @. tends_ρae_tot += - -∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea * aux_up[i].h_tot)) + + -∇c(LBF(Ic(CCG.WVector(w_up)) * ρarea * aux_up.:($$i).h_tot)) + volume_term @. tends_ρaq_tot += -∇c(LBF(Ic(CCG.WVector(w_up)) * ρaq_tot)) @@ -421,8 +422,8 @@ function compute_implicit_up_tendencies!( grad_f = CCO.GradientC2F(; prog_bcs...) @inbounds for i in 1:N_up - w_up = prog_up_f[i].w - tends_w = tendencies_up_f[i].w + w_up = prog_up_f.:($i).w + tends_w = tendencies_up_f.:($i).w @. tends_w += -grad_f(LBC(LA.norm_sqr(CCG.WVector(w_up)) / 2)) tends_w[kf_surf] = zero(tends_w[kf_surf]) end @@ -460,55 +461,58 @@ function compute_explicit_up_tendencies!( @inbounds for i in 1:N_up - w_up = prog_up_f[i].w + w_up = prog_up_f.:($i).w w_en = aux_en_f.w # Augment the tendencies of updraft area, tracers and vertical velocity # entrainment and detrainment - could be moved to implicit volume_term_entr = @. -p_c / ρ_c * - prog_up[i].ρarea * + prog_up.:($$i).ρarea * Ic(wcomponent(CCG.WVector(w_up))) * - (aux_up[i].entr - aux_up[i].detr) - @. tendencies_up[i].ρarea += - prog_up[i].ρarea * + (aux_up.:($$i).entr - aux_up.:($$i).detr) + @. tendencies_up.:($$i).ρarea += + prog_up.:($$i).ρarea * Ic(wcomponent(CCG.WVector(w_up))) * - (aux_up[i].entr - aux_up[i].detr) - @. tendencies_up[i].ρae_tot += - prog_up[i].ρarea * + (aux_up.:($$i).entr - aux_up.:($$i).detr) + @. tendencies_up.:($$i).ρae_tot += + prog_up.:($$i).ρarea * aux_en.h_tot * Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].entr - - prog_up[i].ρarea * - aux_up[i].h_tot * + aux_up.:($$i).entr - + prog_up.:($$i).ρarea * + aux_up.:($$i).h_tot * Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr + volume_term_entr - @. tendencies_up[i].ρaq_tot += - prog_up[i].ρarea * + aux_up.:($$i).detr + volume_term_entr + @. tendencies_up.:($$i).ρaq_tot += + prog_up.:($$i).ρarea * aux_en.q_tot * Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].entr - - prog_up[i].ρaq_tot * + aux_up.:($$i).entr - + prog_up.:($$i).ρaq_tot * Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr - @. tendencies_up_f[i].w += - w_up * I0f(aux_up[i].entr) * (wcomponent(CCG.WVector(w_en - w_up))) + aux_up.:($$i).detr + @. tendencies_up_f.:($$i).w += + w_up * + I0f(aux_up.:($$i).entr) * + (wcomponent(CCG.WVector(w_en - w_up))) # precipitation formation - @. tendencies_up[i].ρae_tot += - prog_gm.ρ * aux_up[i].e_tot_tendency_precip_formation - @. tendencies_up[i].ρaq_tot += - prog_gm.ρ * aux_up[i].qt_tendency_precip_formation + @. tendencies_up.:($$i).ρae_tot += + prog_gm.ρ * aux_up.:($$i).e_tot_tendency_precip_formation + @. tendencies_up.:($$i).ρaq_tot += + prog_gm.ρ * aux_up.:($$i).qt_tendency_precip_formation # buoyancy and pressure - @. tendencies_up_f[i].w += CCG.Covariant3Vector( - CCG.WVector(aux_up_f[i].buoy + aux_up_f[i].nh_pressure), + @. tendencies_up_f.:($$i).w += CCG.Covariant3Vector( + CCG.WVector(aux_up_f.:($$i).buoy + aux_up_f.:($$i).nh_pressure), ) # TODO - to be removed? - tendencies_up[i].ρarea[kc_surf] = 0 - tendencies_up[i].ρae_tot[kc_surf] = 0 - tendencies_up[i].ρaq_tot[kc_surf] = 0 - tendencies_up_f[i].w[kf_surf] = zero(tendencies_up_f[i].w[kf_surf]) + tendencies_up.:($i).ρarea[kc_surf] = 0 + tendencies_up.:($i).ρae_tot[kc_surf] = 0 + tendencies_up.:($i).ρaq_tot[kc_surf] = 0 + tendencies_up_f.:($i).w[kf_surf] = + zero(tendencies_up_f.:($i).w[kf_surf]) end return nothing end @@ -545,22 +549,23 @@ function filter_updraft_vars( a_max = edmf.max_area @inbounds for i in 1:N_up - @. aux_bulk.filter_flag_1 = ifelse(prog_up[i].ρarea < FT(0), 1, 0) - @. aux_bulk.filter_flag_3 = ifelse(prog_up[i].ρaq_tot < FT(0), 1, 0) - @. aux_bulk.filter_flag_4 = ifelse(prog_up[i].ρarea > ρ_c * a_max, 1, 0) - - @. prog_up[i].ρarea = max(prog_up[i].ρarea, 0) #flag_1 - @. prog_up[i].ρaq_tot = max(prog_up[i].ρaq_tot, 0) #flag_3 - @. prog_up[i].ρarea = min(prog_up[i].ρarea, ρ_c * a_max) #flag_4 + @. aux_bulk.filter_flag_1 = ifelse(prog_up.:($$i).ρarea < FT(0), 1, 0) + @. aux_bulk.filter_flag_3 = ifelse(prog_up.:($$i).ρaq_tot < FT(0), 1, 0) + @. aux_bulk.filter_flag_4 = + ifelse(prog_up.:($$i).ρarea > ρ_c * a_max, 1, 0) + + @. prog_up.:($$i).ρarea = max(prog_up.:($$i).ρarea, 0) #flag_1 + @. prog_up.:($$i).ρaq_tot = max(prog_up.:($$i).ρaq_tot, 0) #flag_3 + @. prog_up.:($$i).ρarea = min(prog_up.:($$i).ρarea, ρ_c * a_max) #flag_4 end @inbounds for i in 1:N_up - @. prog_up_f[i].w = CCG.Covariant3Vector( - CCG.WVector(max(wcomponent(CCG.WVector(prog_up_f[i].w)), 0)), + @. prog_up_f.:($$i).w = CCG.Covariant3Vector( + CCG.WVector(max(wcomponent(CCG.WVector(prog_up_f.:($$i).w)), 0)), ) a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) If = CCO.InterpolateC2F(; a_up_bcs...) - @. prog_up_f[i].w = - Int(If(prog_up[i].ρarea) >= ρ_f * a_min) * prog_up_f[i].w + @. prog_up_f.:($$i).w = + Int(If(prog_up.:($$i).ρarea) >= ρ_f * a_min) * prog_up_f.:($$i).w end Δz = Fields.Δz_field(axes(ρ_c)) @@ -570,51 +575,59 @@ function filter_updraft_vars( # this is needed to make sure Rico is unchanged. # TODO : look into it further to see why # a similar filtering of ρaθ_liq_ice breaks the simulation - @. prog_up[i].ρaq_tot = ifelse( + @. prog_up.:($$i).ρaq_tot = ifelse( z > Δz1 / 2, ifelse( - prog_up[i].ρarea / ρ_c < a_min, + prog_up.:($$i).ρarea / ρ_c < a_min, 0, - max(prog_up[i].ρaq_tot, 0), + max(prog_up.:($$i).ρaq_tot, 0), ), - prog_up[i].ρaq_tot, + prog_up.:($$i).ρaq_tot, ) - @. prog_up[i].ρarea = ifelse( + @. prog_up.:($$i).ρarea = ifelse( z > Δz1 / 2, - ifelse(prog_up[i].ρarea / ρ_c < a_min, 0, max(prog_up[i].ρarea, 0)), - prog_up[i].ρarea, + ifelse( + prog_up.:($$i).ρarea / ρ_c < a_min, + 0, + max(prog_up.:($$i).ρarea, 0), + ), + prog_up.:($$i).ρarea, ) - @. prog_up[i].ρae_tot = ifelse( + @. prog_up.:($$i).ρae_tot = ifelse( z > Δz1 / 2, - ifelse(prog_up[i].ρarea / ρ_c < a_min, 0, prog_up[i].ρae_tot), - prog_up[i].ρae_tot, + ifelse( + prog_up.:($$i).ρarea / ρ_c < a_min, + 0, + prog_up.:($$i).ρae_tot, + ), + prog_up.:($$i).ρae_tot, ) end Ic = CCO.InterpolateF2C() C123 = CCG.Covariant123Vector @inbounds for i in 1:N_up - @. prog_up[i].ρarea = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, + @. prog_up.:($$i).ρarea = ifelse( + Ic(wcomponent(CCG.WVector(prog_up_f.:($$i).w))) <= 0, FT(0), - prog_up[i].ρarea, + prog_up.:($$i).ρarea, ) - @. prog_up[i].ρae_tot = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, + @. prog_up.:($$i).ρae_tot = ifelse( + Ic(wcomponent(CCG.WVector(prog_up_f.:($$i).w))) <= 0, FT(0), - prog_up[i].ρae_tot, + prog_up.:($$i).ρae_tot, ) - @. prog_up[i].ρaq_tot = ifelse( - Ic(wcomponent(CCG.WVector(prog_up_f[i].w))) <= 0, + @. prog_up.:($$i).ρaq_tot = ifelse( + Ic(wcomponent(CCG.WVector(prog_up_f.:($$i).w))) <= 0, FT(0), - prog_up[i].ρaq_tot, + prog_up.:($$i).ρaq_tot, ) θ_surf = θ_surface_bc(grid, state, edmf, i, param_set) q_surf = q_surface_bc(grid, state, edmf, i, param_set) a_surf = area_surface_bc(surface_conditions, edmf) e_kin = @. LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), + C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f.:($$i).w))), ) / 2 e_pot_surf = geopotential(thermo_params, grid.zc.z[kc_surf]) ts_up_i_surf = @@ -625,9 +638,10 @@ function filter_updraft_vars( e_kin[kc_surf], e_pot_surf, ) - prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf - prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf + prog_up.:($i).ρarea[kc_surf] = ρ_c[kc_surf] * a_surf + prog_up.:($i).ρae_tot[kc_surf] = + prog_up.:($i).ρarea[kc_surf] * e_tot_surf + prog_up.:($i).ρaq_tot[kc_surf] = prog_up.:($i).ρarea[kc_surf] * q_surf end return nothing end diff --git a/src/TurbulenceConvection_deprecated/Fields.jl b/src/TurbulenceConvection_deprecated/Fields.jl index 9319e5c523e..cbcef2d4414 100644 --- a/src/TurbulenceConvection_deprecated/Fields.jl +++ b/src/TurbulenceConvection_deprecated/Fields.jl @@ -34,35 +34,16 @@ const CenterFields = Union{ CC.Fields.CenterFiniteDifferenceField, } -Base.@propagate_inbounds Base.getindex(field::FDFields, i::Integer) = - Base.getproperty(field, i) - Base.@propagate_inbounds Base.getindex(field::CenterFields, i::Cent) = Base.getindex(CC.Fields.field_values(field), i.i) Base.@propagate_inbounds Base.setindex!(field::CenterFields, v, i::Cent) = Base.setindex!(CC.Fields.field_values(field), v, i.i) -Base.@propagate_inbounds Base.getindex(field::FaceFields, i::CCO.PlusHalf) = - Base.getindex(CC.Fields.field_values(field), i.i) -Base.@propagate_inbounds Base.setindex!(field::FaceFields, v, i::CCO.PlusHalf) = - Base.setindex!(CC.Fields.field_values(field), v, i.i) - Base.@propagate_inbounds Base.getindex(field::FaceFields, ::Cent) = error("Attempting to getindex with a center index (Cent) into a Face field") -Base.@propagate_inbounds Base.getindex(field::CenterFields, ::CCO.PlusHalf) = - error( - "Attempting to getindex with a face index (PlusHalf) into a Center field", - ) Base.@propagate_inbounds Base.setindex!(field::FaceFields, v, ::Cent) = error("Attempting to setindex with a center index (Cent) into a Face field") -Base.@propagate_inbounds Base.setindex!( - field::CenterFields, - v, - ::CCO.PlusHalf, -) = error( - "Attempting to setindex with a face index (PlusHalf) into a Center field", -) function Base.cumsum(field::CC.Fields.FiniteDifferenceField) Base.cumsum( diff --git a/src/TurbulenceConvection_deprecated/update_aux.jl b/src/TurbulenceConvection_deprecated/update_aux.jl index 91ecab92767..647bdeb01ae 100644 --- a/src/TurbulenceConvection_deprecated/update_aux.jl +++ b/src/TurbulenceConvection_deprecated/update_aux.jl @@ -67,56 +67,64 @@ function update_aux!( C123 = CCG.Covariant123Vector @inbounds for i in 1:N_up - @. aux_up[i].e_kin = + @. aux_up.:($$i).e_kin = LA.norm_sqr( - C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f[i].w))), + C123(prog_gm_uₕ) + C123(Ic(CCG.WVector(prog_up_f.:($$i).w))), ) / 2 end e_pot(z) = geopotential(thermo_params, z) - thresh_area(prog_up, ρ_c) = prog_up[i].ρarea / ρ_c[k] >= edmf.minimum_area + @inbounds for i in 1:N_up - @. aux_up[i].e_tot = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρae_tot / prog_up[i].ρarea, + @. aux_up.:($$i).e_tot = ifelse( + prog_up.:($$i).ρarea / ρ_c >= edmf.minimum_area, + prog_up.:($$i).ρae_tot / prog_up.:($$i).ρarea, aux_gm.e_tot, ) - @. aux_up[i].q_tot = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρaq_tot / prog_up[i].ρarea, + @. aux_up.:($$i).q_tot = ifelse( + prog_up.:($$i).ρarea / ρ_c >= edmf.minimum_area, + prog_up.:($$i).ρaq_tot / prog_up.:($$i).ρarea, aux_gm.q_tot, ) - @. aux_up[i].area = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - prog_up[i].ρarea / ρ_c, + @. aux_up.:($$i).area = ifelse( + prog_up.:($$i).ρarea / ρ_c >= edmf.minimum_area, + prog_up.:($$i).ρarea / ρ_c, 0, ) - @. aux_up[i].e_kin = ifelse( - prog_up[i].ρarea / ρ_c >= edmf.minimum_area, - aux_up[i].e_kin, + @. aux_up.:($$i).e_kin = ifelse( + prog_up.:($$i).ρarea / ρ_c >= edmf.minimum_area, + aux_up.:($$i).e_kin, e_kin, ) ##### ##### Set primitive variables ##### - e_int = @. aux_up[i].e_tot - aux_up[i].e_kin - e_pot(zc) + e_int = @. aux_up.:($$i).e_tot - aux_up.:($$i).e_kin - e_pot(zc) if edmf.moisture_model isa DryModel - @. aux_up[i].ts = TD.PhaseDry_pe(thermo_params, p_c, e_int) + @. aux_up.:($$i).ts = TD.PhaseDry_pe(thermo_params, p_c, e_int) elseif edmf.moisture_model isa EquilMoistModel - @. aux_up[i].ts = - TD.PhaseEquil_peq(thermo_params, p_c, e_int, aux_up[i].q_tot) + @. aux_up.:($$i).ts = TD.PhaseEquil_peq( + thermo_params, + p_c, + e_int, + aux_up.:($$i).q_tot, + ) elseif edmf.moisture_model isa NonEquilMoistModel error("Unsupported moisture model") end - ts_up = aux_up[i].ts - @. aux_up[i].θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts_up) - @. aux_up[i].h_tot = - TD.total_specific_enthalpy(thermo_params, ts_up, aux_up[i].e_tot) - @. aux_up[i].q_liq = TD.liquid_specific_humidity(thermo_params, ts_up) - @. aux_up[i].q_ice = TD.ice_specific_humidity(thermo_params, ts_up) - @. aux_up[i].T = TD.air_temperature(thermo_params, ts_up) - @. aux_up[i].RH = TD.relative_humidity(thermo_params, ts_up) + ts_up = aux_up.:($i).ts + @. aux_up.:($$i).θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ts_up) + @. aux_up.:($$i).h_tot = TD.total_specific_enthalpy( + thermo_params, + ts_up, + aux_up.:($$i).e_tot, + ) + @. aux_up.:($$i).q_liq = + TD.liquid_specific_humidity(thermo_params, ts_up) + @. aux_up.:($$i).q_ice = TD.ice_specific_humidity(thermo_params, ts_up) + @. aux_up.:($$i).T = TD.air_temperature(thermo_params, ts_up) + @. aux_up.:($$i).RH = TD.relative_humidity(thermo_params, ts_up) end ##### @@ -124,14 +132,16 @@ function update_aux!( ##### @. aux_bulk.area = 0 @inbounds for i in 1:N_up - @. aux_bulk.area += aux_up[i].area + @. aux_bulk.area += aux_up.:($$i).area end @. aux_bulk.q_tot = 0 @. aux_bulk.h_tot = 0 @inbounds for i in 1:N_up - @. aux_bulk.q_tot += aux_up[i].area * aux_up[i].q_tot / aux_bulk.area - @. aux_bulk.h_tot += aux_up[i].area * aux_up[i].h_tot / aux_bulk.area + @. aux_bulk.q_tot += + aux_up.:($$i).area * aux_up.:($$i).q_tot / aux_bulk.area + @. aux_bulk.h_tot += + aux_up.:($$i).area * aux_up.:($$i).h_tot / aux_bulk.area end @inbounds for i in 1:N_up @. aux_bulk.q_tot = @@ -146,7 +156,8 @@ function update_aux!( @. aux_en_f.w = prog_gm_f.u₃ / (1 - Ifb(aux_bulk.area)) @inbounds for i in 1:N_up @. aux_en_f.w -= - Ifb(aux_up[i].area) * prog_up_f[i].w / (1 - Ifb(aux_bulk.area)) + Ifb(aux_up.:($$i).area) * prog_up_f.:($$i).w / + (1 - Ifb(aux_bulk.area)) end @. aux_en.e_kin = @@ -265,18 +276,18 @@ function update_aux!( LBF_ρ(TD.air_density(thermo_params, aux_en.ts)), ) @inbounds for i in 1:N_up - @. aux_up_f[i].buoy = buoyancy_c( + @. aux_up_f.:($$i).buoy = buoyancy_c( thermo_params, ρ_f, - LBF_ρ(TD.air_density(thermo_params, aux_up[i].ts)), + LBF_ρ(TD.air_density(thermo_params, aux_up.:($$i).ts)), ) end @. aux_gm_f.buoy = (1.0 - LBF_a(aux_bulk.area)) * aux_en_f.buoy @inbounds for i in 1:N_up - @. aux_gm_f.buoy += LBF_a(aux_up[i].area) * aux_up_f[i].buoy + @. aux_gm_f.buoy += LBF_a(aux_up.:($$i).area) * aux_up_f.:($$i).buoy end @inbounds for i in 1:N_up - @. aux_up_f[i].buoy -= aux_gm_f.buoy + @. aux_up_f.:($$i).buoy -= aux_gm_f.buoy end # Only needed for diagnostics/plotting @. aux_en_f.buoy -= aux_gm_f.buoy @@ -290,9 +301,11 @@ function update_aux!( @. aux_bulk.q_ice = 0 @. aux_bulk.T = 0 @inbounds for i in 1:N_up - @. aux_bulk.q_liq += aux_up[i].area * aux_up[i].q_liq / aux_bulk.area - @. aux_bulk.q_ice += aux_up[i].area * aux_up[i].q_ice / aux_bulk.area - @. aux_bulk.T += aux_up[i].area * aux_up[i].T / aux_bulk.area + @. aux_bulk.q_liq += + aux_up.:($$i).area * aux_up.:($$i).q_liq / aux_bulk.area + @. aux_bulk.q_ice += + aux_up.:($$i).area * aux_up.:($$i).q_ice / aux_bulk.area + @. aux_bulk.T += aux_up.:($$i).area * aux_up.:($$i).T / aux_bulk.area end @. aux_bulk.q_liq = ifelse(aux_bulk.area > 0, aux_bulk.q_liq, 0) @. aux_bulk.q_ice = ifelse(aux_bulk.area > 0, aux_bulk.q_ice, 0) @@ -322,8 +335,8 @@ function update_aux!( @inbounds for i in 1:N_up @. aux_gm.tke += 0.5 * - aux_up[i].area * - LA.norm_sqr(C123(Ic(wvec(prog_up_f[i].w - prog_gm_f.u₃)))) + aux_up.:($$i).area * + LA.norm_sqr(C123(Ic(wvec(prog_up_f.:($$i).w - prog_gm_f.u₃)))) end ##### @@ -331,12 +344,12 @@ function update_aux!( ##### parent(aux_tc_f.bulk.w) .= 0 @inbounds for i in 1:N_up - a_up = aux_up[i].area + a_up = aux_up.:($i).area a_up_bcs = a_up_boundary_conditions(surface_conditions, edmf) Ifu = CCO.InterpolateC2F(; a_up_bcs...) @. aux_tc_f.bulk.w += ifelse( Ifb(aux_bulk.area) > 0, - Ifu(a_up) * prog_up_f[i].w / Ifb(aux_bulk.area), + Ifu(a_up) * prog_up_f.:($$i).w / Ifb(aux_bulk.area), CCG.Covariant3Vector(FT(0)), ) end @@ -346,11 +359,11 @@ function update_aux!( ##### # entrainment and detrainment @inbounds for i in 1:N_up - @. aux_up[i].entr = FT(5e-4) - @. aux_up[i].detr = pi_groups_detrainment!( + @. aux_up.:($$i).entr = FT(5e-4) + @. aux_up.:($$i).detr = pi_groups_detrainment!( aux_gm.tke, - aux_up[i].area, - aux_up[i].RH, + aux_up.:($$i).area, + aux_up.:($$i).RH, aux_en.area, aux_en.tke, aux_en.RH, @@ -361,13 +374,13 @@ function update_aux!( (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0)))) ∇ = CCO.DivergenceC2F(; w_bcs...) @inbounds for i in 1:N_up - updraft_top = compute_updraft_top(aux_up[i].area) - @. aux_up_f[i].nh_pressure = compute_nh_pressure!( + updraft_top = compute_updraft_top(aux_up.:($i).area) + @. aux_up_f.:($$i).nh_pressure = compute_nh_pressure!( param_set, - aux_up_f[i].buoy, - wcomponent(CCG.WVector(prog_up_f[i].w)), - ∇(Ic(CCG.WVector(prog_up_f[i].w))), - wcomponent(CCG.WVector(prog_up_f[i].w - aux_en_f.w)), + aux_up_f.:($$i).buoy, + wcomponent(CCG.WVector(prog_up_f.:($$i).w)), + ∇(Ic(CCG.WVector(prog_up_f.:($$i).w))), + wcomponent(CCG.WVector(prog_up_f.:($$i).w - aux_en_f.w)), updraft_top, ) end @@ -382,11 +395,11 @@ function update_aux!( parent(b_exch) .= 0 w_en = aux_en_f.w @inbounds for i in 1:N_up - w_up = prog_up_f[i].w + w_up = prog_up_f.:($i).w @. b_exch += - aux_up[i].area * + aux_up.:($$i).area * Ic(wcomponent(CCG.WVector(w_up))) * - aux_up[i].detr / aux_en.area * + aux_up.:($$i).detr / aux_en.area * (1 / 2 * (Ic(wcomponent(CCG.WVector(w_up - w_en))))^2 - aux_en.tke) end @@ -526,31 +539,31 @@ function update_aux!( param_set, edmf.precip_model, prog_gm, - aux_up[i].area, + aux_up.:($$i).area, zc, Δt, - aux_up[i].ts, + aux_up.:($$i).ts, ) - @. aux_up[i].θ_liq_ice_tendency_precip_formation = - mph.θ_liq_ice_tendency * aux_up[i].area - @. aux_up[i].e_tot_tendency_precip_formation = - mph.e_tot_tendency * aux_up[i].area - @. aux_up[i].qt_tendency_precip_formation = - mph.qt_tendency * aux_up[i].area - @. aux_up[i].qr_tendency_precip_formation = - mph.qr_tendency * aux_up[i].area - @. aux_up[i].qs_tendency_precip_formation = - mph.qs_tendency * aux_up[i].area + @. aux_up.:($$i).θ_liq_ice_tendency_precip_formation = + mph.θ_liq_ice_tendency * aux_up.:($$i).area + @. aux_up.:($$i).e_tot_tendency_precip_formation = + mph.e_tot_tendency * aux_up.:($$i).area + @. aux_up.:($$i).qt_tendency_precip_formation = + mph.qt_tendency * aux_up.:($$i).area + @. aux_up.:($$i).qr_tendency_precip_formation = + mph.qr_tendency * aux_up.:($$i).area + @. aux_up.:($$i).qs_tendency_precip_formation = + mph.qs_tendency * aux_up.:($$i).area # TODO - to be deleted once we sum all tendencies elsewhere @. aux_bulk.e_tot_tendency_precip_formation += - aux_up[i].e_tot_tendency_precip_formation + aux_up.:($$i).e_tot_tendency_precip_formation @. aux_bulk.qt_tendency_precip_formation += - aux_up[i].qt_tendency_precip_formation + aux_up.:($$i).qt_tendency_precip_formation @. aux_bulk.qr_tendency_precip_formation += - aux_up[i].qr_tendency_precip_formation + aux_up.:($$i).qr_tendency_precip_formation @. aux_bulk.qs_tendency_precip_formation += - aux_up[i].qs_tendency_precip_formation + aux_up.:($$i).qs_tendency_precip_formation end return nothing diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 88e8e95d0cd..2be28412bcd 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -358,9 +358,9 @@ values of the first updraft. function output_diagnostic_sgs_quantities(Y, p, t) thermo_params = CAP.thermodynamics_params(p.params) (; ᶜρaʲs, ᶜtsʲs) = p - ᶠu³⁺ = p.ᶠu³ʲs[1] + ᶠu³⁺ = p.ᶠu³ʲs.:1 ᶜu⁺ = @. (C123(Y.c.uₕ) + C123(ᶜinterp(ᶠu³⁺))) - ᶜts⁺ = @. ᶜtsʲs[1] - ᶜa⁺ = @. draft_area(ᶜρaʲs[1], TD.air_density(thermo_params, ᶜts⁺)) + ᶜts⁺ = @. ᶜtsʲs.:1 + ᶜa⁺ = @. draft_area(ᶜρaʲs.:1, TD.air_density(thermo_params, ᶜts⁺)) return (; ᶜu⁺, ᶠu³⁺, ᶜts⁺, ᶜa⁺) end