diff --git a/src/cache/advective_edmf_precomputed_quantities.jl b/src/cache/advective_edmf_precomputed_quantities.jl index 96067490247..e87bf1d37af 100644 --- a/src/cache/advective_edmf_precomputed_quantities.jl +++ b/src/cache/advective_edmf_precomputed_quantities.jl @@ -25,15 +25,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p - (; - ᶜuʲs, - ᶠu³ʲs, - ᶜKʲs, - ᶜtsʲs, - ᶜρʲs, - ᶜentrʲs, - ᶜdetrʲs, - ) = p + (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions @. ᶜρa⁰ = ρa⁰(Y.c) @@ -56,12 +48,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model) set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³) @. ᶜK⁰ += ᶜtke⁰ - @. ᶜts⁰ = TD.PhaseEquil_phq( - thermo_params, - ᶜp, - ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ, - ᶜq_tot⁰, - ) + @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ, ᶜq_tot⁰) @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰) for j in 1:n @@ -75,13 +62,8 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t) ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³) - # FIXME: get ts - @. ᶜtsʲ = TD.PhaseEquil_phq( - thermo_params, - ᶜp, - ᶜh_totʲ - ᶜKʲ - ᶜΦ, - ᶜq_totʲ, - ) + @. ᶜtsʲ = + TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_totʲ - ᶜKʲ - ᶜΦ, ᶜq_totʲ) @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ) # When ᶜe_intʲ = ᶜe_int and ᶜq_totʲ = ᶜq_tot, we still observe that diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index f22c36482e0..e0b34a1a12f 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -147,7 +147,12 @@ function precomputed_quantities(Y, atmos) ᶜK_h = similar(Y.c, FT), ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}), ) : (;) - return (; gs_quantities..., sgs_quantities..., advective_sgs_quantities..., diagnostic_sgs_quantities...) + return (; + gs_quantities..., + sgs_quantities..., + advective_sgs_quantities..., + diagnostic_sgs_quantities..., + ) end # Interpolates the third contravariant component of Y.c.uₕ to cell faces. diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl index 2ad588ef346..8caf9e21705 100644 --- a/src/initial_conditions/atmos_state.jl +++ b/src/initial_conditions/atmos_state.jl @@ -144,9 +144,9 @@ function turbconv_center_variables(ls, turbconv_model::AdvectiveEDMFX, gs_vars) sgs⁰ = (; ρatke = ls.ρ * (1 - a_draft) * ls.turbconv_state.tke) ρa = ls.ρ * a_draft / n h_tot = - TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) + - norm_sqr(ls.velocity) / 2 + - CAP.grav(ls.params) * ls.geometry.coordinates.z + TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) + + norm_sqr(ls.velocity) / 2 + + CAP.grav(ls.params) * ls.geometry.coordinates.z q_tot = TD.total_specific_humidity(ls.thermo_params, ls.thermo_state) sgsʲs = ntuple(_ -> (; ρa = ρa, h_tot = h_tot, q_tot = q_tot), Val(n)) return (; sgs⁰, sgsʲs) diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 0bc926cdcff..dcfe9979b5f 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -12,12 +12,14 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) if p.atmos.turbconv_model isa AbstractEDMF (; ᶜu⁰) = p end - if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX (; ᶜuʲs) = p end @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu) - if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX for j in 1:n @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j)) end @@ -108,7 +110,9 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) (; ᶠu³⁰) = advect_tke ? p : all_nothing ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing - ᶜtke⁰ = advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing + ᶜtke⁰ = + advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : + nothing ᶜa_scalar = p.ᶜtemp_scalar ᶜω³ = p.ᶜtemp_CT3 ᶠω¹² = p.ᶠtemp_CT12 @@ -173,8 +177,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j)) @. ᶜa_scalar[colidx] = - ᶜh_totʲs.:($$j)[colidx] * - draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx]) + ᶜh_totʲs.:($$j)[colidx] * draft_area( + Y.c.sgsʲs.:($$j).ρa[colidx], + ᶜρʲs.:($$j)[colidx], + ) vertical_transport!( Yₜ.c.sgsʲs.:($j).ρae_tot[colidx], ᶜJ[colidx], @@ -190,8 +196,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j)) χ_name == :e_tot && continue @. ᶜa_scalar[colidx] = - ᶜχʲ[colidx] * - draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx]) + ᶜχʲ[colidx] * draft_area( + Y.c.sgsʲs.:($$j).ρa[colidx], + ᶜρʲs.:($$j)[colidx], + ) vertical_transport!( ᶜρaχʲₜ[colidx], ᶜJ[colidx], @@ -220,11 +228,15 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) edmfx_upwinding, ) - if :h_tot in propertynames(Yₜ.c.sgsʲs.:($j)) - @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] += 0 - end - - @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] += 0 + # FIXME: boundary conditions + @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] -= ᶜinterp( + adjoint(ᶠu³ʲs.:($$j)[colidx]) * + ᶠgradᵥ(Y.c.sgsʲs.:($$j).h_tot[colidx]), + ) + @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] -= ᶜinterp( + adjoint(ᶠu³ʲs.:($$j)[colidx]) * + ᶠgradᵥ(Y.c.sgsʲs.:($$j).q_tot[colidx]), + ) end end diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index 27e2982ffce..90b27fe0fde 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -315,7 +315,14 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX) return nothing end -function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX) +function edmfx_entr_detr_tendency!( + Yₜ, + Y, + p, + t, + colidx, + turbconv_model::AdvectiveEDMFX, +) n = n_mass_flux_subdomains(turbconv_model) (; ᶜentrʲs, ᶜdetrʲs) = p @@ -330,7 +337,7 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::Advect @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] += ᶜentrʲs.:($$j)[colidx] * (ᶜh_tot⁰[colidx] - Y.c.sgsʲs.:($$j).h_tot[colidx]) - + @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] += ᶜentrʲs.:($$j)[colidx] * (ᶜq_tot⁰[colidx] - Y.c.sgsʲs.:($$j).q_tot[colidx]) diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl index 427964387c2..5a3f5189ac6 100644 --- a/src/prognostic_equations/edmfx_sgs_flux.jl +++ b/src/prognostic_equations/edmfx_sgs_flux.jl @@ -90,8 +90,14 @@ function edmfx_sgs_mass_flux_tendency!( end # FIXME: add sgs mass flux -edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX) = - nothing +edmfx_sgs_mass_flux_tendency!( + Yₜ, + Y, + p, + t, + colidx, + turbconv_model::AdvectiveEDMFX, +) = nothing function edmfx_sgs_mass_flux_tendency!( Yₜ, @@ -225,7 +231,14 @@ function edmfx_sgs_diffusive_flux_tendency!( end # FIXME: add SGS diffusive flux -edmfx_sgs_diffusive_flux_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX) = nothing +edmfx_sgs_diffusive_flux_tendency!( + Yₜ, + Y, + p, + t, + colidx, + turbconv_model::AdvectiveEDMFX, +) = nothing function edmfx_sgs_diffusive_flux_tendency!( Yₜ, diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index dbb1d19105b..a323383d09d 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -178,7 +178,7 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) end end # FIXME: add hyperdiffusion for advective edmf - + if turbconv_model isa DiagnosticEDMFX && diffuse_tke @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰)) end diff --git a/src/prognostic_equations/zero_velocity.jl b/src/prognostic_equations/zero_velocity.jl index 4f926c26f5a..cb6787ddeb1 100644 --- a/src/prognostic_equations/zero_velocity.jl +++ b/src/prognostic_equations/zero_velocity.jl @@ -10,7 +10,8 @@ function zero_velocity_tendency!(Yₜ, Y, p, t, colidx) @. Yₜ.c.uₕ[colidx] = C12(FT(0), FT(0)) @. Yₜ.f.u₃[colidx] = Geometry.Covariant3Vector(FT(0)) - if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX + if p.atmos.turbconv_model isa EDMFX || + p.atmos.turbconv_model isa AdvectiveEDMFX for j in 1:n @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] = Geometry.Covariant3Vector(FT(0)) diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index ecdff9007cb..16df9994182 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -341,7 +341,8 @@ function get_turbconv_model( turbconv_params, ) turbconv = parsed_args["turbconv"] - @assert turbconv in (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx") + @assert turbconv in + (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx") return if turbconv == "edmf" TC.EDMFModel( diff --git a/src/solver/types.jl b/src/solver/types.jl index 92ed541eda6..26f58441d39 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -160,7 +160,8 @@ EDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = EDMFX{N, TKE, FT}(a_half) struct AdvectiveEDMFX{N, TKE, FT} <: AbstractEDMF a_half::FT # WARNING: this should never be used outside of divide_by_ρa end -AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = AdvectiveEDMFX{N, TKE, FT}(a_half) +AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = + AdvectiveEDMFX{N, TKE, FT}(a_half) struct DiagnosticEDMFX{N, TKE, FT} <: AbstractEDMF a_int::FT # area fraction of the first interior cell above the surface diff --git a/toml/edmfx_box.toml b/toml/edmfx_box.toml index 4f58178f41d..5bdb7ee68ab 100644 --- a/toml/edmfx_box.toml +++ b/toml/edmfx_box.toml @@ -13,3 +13,8 @@ alias = "min_area" value = 1.0e-3 type = "float" description = "Minimum area fraction per updraft. Parameter not described in the literature." + +[entr_coeff] +alias = "entr_coeff" +value = 1.0e-1 +type = "float" \ No newline at end of file