remove theta as an energy variable

[szy21]

Purpose

bye bye theta and energy form

To-do

Content

---

• I have read and checked the items on the review checklist.

[szy21]

I need to remove the ρθ option in the implicit solver too. I'll do it now.

[Sbozzolo]

• The default config should be updated to remove the old option.
• If you remove energy_form from the AtmosModel, all the references to TotalEnergy should disappear. There's a handful left (most importantly, in `default_diagnostics)
• Some of the post processing steps have references to the energy variable
• There's still many ρθs around

❯ grep -r "ρθ"
post_processing/remap/remap_helpers.jl:    elseif :ρθ in propertynames(Y.c)
post_processing/remap/remap_helpers.jl:    elseif :ρθ in propertynames(Y.c)
post_processing/remap/remap_helpers.jl:        nc_thermo[:, 1] = Y.c.ρθ ./ Y.c.ρ
src/cache/precomputed_quantities.jl:        TD.PhaseDry_ρθ(thermo_params, ρ, θ)
src/cache/precomputed_quantities.jl:        TD.PhaseEquil_ρθq(thermo_params, ρ, θ, q_tot)
src/cache/precomputed_quantities.jl:        TD.PhaseNonEquil_ρθq(thermo_params, ρ, θ, q_pt)
src/parameterized_tendencies/radiation/held_suarez.jl:    if :ρθ in propertynames(Y.c)
src/parameterized_tendencies/radiation/held_suarez.jl:        @. Yₜ.c.ρθ[colidx] -=
src/parameterized_tendencies/radiation/radiation.jl:    if :ρθ in propertynames(Y.c)
src/parameterized_tendencies/radiation/radiation.jl:        error("radiation_tendency! not implemented for ρθ")
src/prognostic_equations/implicit/implicit_solver.jl:        MatrixFields.unrolled_findonly(is_in_Y, (@name(c.ρe_tot), @name(c.ρθ)))
src/prognostic_equations/implicit/implicit_solver.jl:    # ᶜp = p_ref_theta * (ᶜρθ * R_d / p_ref_theta)^(1 / (1 - κ_d)) or
src/prognostic_equations/implicit/implicit_solver.jl:    # ∂(ᶜp)/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #         R_d / (1 - κ_d) * (ᶜρθ * R_d / p_ref_theta)^(κ_d / (1 - κ_d))
src/prognostic_equations/implicit/implicit_solver.jl:        (@name(c.ρθ), @name(ᶜspecific.θ), energy_upwinding),
src/prognostic_equations/implicit/implicit_solver.jl:    # ∂(ᶠu₃ₜ)/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #     ∂(ᶠu₃ₜ)/∂(ᶠgradᵥ(ᶜp - ᶜp_ref)) ⋅ ∂(ᶠgradᵥ(ᶜp - ᶜp_ref))/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #     ∂(ᶜp)/∂(ᶜρθ) and
src/prognostic_equations/implicit/implicit_solver.jl:    if MatrixFields.has_field(Y, @name(c.ρθ))
src/prognostic_equations/implicit/implicit_solver.jl:        ᶜρθ = Y.c.ρθ
src/prognostic_equations/implicit/implicit_solver.jl:        ∂ᶠu₃_err_∂ᶜρθ = matrix[@name(f.u₃), @name(c.ρθ)]
src/prognostic_equations/implicit/implicit_solver.jl:        @. ∂ᶠu₃_err_∂ᶜρθ[colidx] =
src/prognostic_equations/implicit/implicit_solver.jl:                (ᶜρθ[colidx] * R_d / p_ref_theta)^(κ_d / (1 - κ_d)),
src/prognostic_equations/advection.jl:    if :ρθ in propertynames(Y.c)
src/prognostic_equations/advection.jl:        @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu)
src/prognostic_equations/hyperdiffusion.jl:    ᶜρ_energyₜ = :θ in propertynames(ᶜspecific) ? Yₜ.c.ρθ : Yₜ.c.ρe_tot
src/utils/utilities.jl:is_energy_var(symbol) = symbol in (:ρθ, :ρe_tot, :ρaθ, :ρae_tot)

Also some :θ

❯ grep -r ":θ"
src/prognostic_equations/implicit/implicit_tendency.jl:            χ_name == :θ ? energy_upwinding : tracer_upwinding,
src/prognostic_equations/hyperdiffusion.jl:    if :θ in propertynames(ᶜspecific)
src/prognostic_equations/hyperdiffusion.jl:    ᶜρ_energyₜ = :θ in propertynames(ᶜspecific) ? Yₜ.c.ρθ : Yₜ.c.ρe_tot
src/prognostic_equations/vertical_diffusion_boundary_layer.jl:        elseif χ_name == :θ
src/utils/variable_manipulations.jl:    Base.structdiff(specific_state, NamedTuple{(:θ, :e_tot)})

In general (there might be some false positives)

❯ grep -rI "θ"
post_processing/plot/plot_helpers.jl:            global θ = Array(nc["potential_temperature"])
post_processing/plot/plot_helpers.jl:            θ = cat(θ, Array(nc["potential_temperature"]), dims = 4)
post_processing/plot/plot_helpers.jl:    θ_timeave_zonalave = calc_zonalave_timeave(θ)
post_processing/plot/plot_helpers.jl:        Z = θ_timeave_zonalave,
post_processing/plot/plot_helpers.jl:        title = "Potential Temperature θ",
post_processing/remap/remap_helpers.jl:    elseif :ρθ in propertynames(Y.c)
post_processing/remap/remap_helpers.jl:    nc_θ = defVar(nc, "potential_temperature", FT, cspace, ("time",))
post_processing/remap/remap_helpers.jl:    elseif :ρθ in propertynames(Y.c)
post_processing/remap/remap_helpers.jl:        nc_thermo[:, 1] = Y.c.ρθ ./ Y.c.ρ
post_processing/remap/remap_helpers.jl:    nc_θ[:, 1] = diag.potential_temperature
src/cache/diagnostic_edmf_precomputed_quantities.jl:            TD.dry_pottemp(thermo_params, ᶜts),                               # θ_sat
src/cache/diagnostic_edmf_precomputed_quantities.jl:            TD.liquid_ice_pottemp(thermo_params, ᶜts),                        # θ_liq_ice_sat
src/cache/diagnostic_edmf_precomputed_quantities.jl:            ),                                                                 # ∂θv∂z_unsat
src/cache/diagnostic_edmf_precomputed_quantities.jl:            ),                                                                 # ∂θl∂z_sat
src/cache/precomputed_quantities.jl:    θ = nothing,
src/cache/precomputed_quantities.jl:    get_ts(ρ::Real, ::Nothing, θ::Real, ::Nothing, ::Nothing, ::Nothing) =
src/cache/precomputed_quantities.jl:        TD.PhaseDry_ρθ(thermo_params, ρ, θ)
src/cache/precomputed_quantities.jl:    get_ts(ρ::Real, ::Nothing, θ::Real, ::Nothing, q_tot::Real, ::Nothing) =
src/cache/precomputed_quantities.jl:        TD.PhaseEquil_ρθq(thermo_params, ρ, θ, q_tot)
src/cache/precomputed_quantities.jl:    get_ts(ρ::Real, ::Nothing, θ::Real, ::Nothing, ::Nothing, q_pt) =
src/cache/precomputed_quantities.jl:        TD.PhaseNonEquil_ρθq(thermo_params, ρ, θ, q_pt)
src/cache/precomputed_quantities.jl:    get_ts(::Nothing, p::Real, θ::Real, ::Nothing, ::Nothing, ::Nothing) =
src/cache/precomputed_quantities.jl:        TD.PhaseDry_pθ(thermo_params, p, θ)
src/cache/precomputed_quantities.jl:    get_ts(::Nothing, p::Real, θ::Real, ::Nothing, q_tot::Real, ::Nothing) =
src/cache/precomputed_quantities.jl:        TD.PhaseEquil_pθq(thermo_params, p, θ, q_tot)
src/cache/precomputed_quantities.jl:    get_ts(::Nothing, p::Real, θ::Real, ::Nothing, ::Nothing, q_pt) =
src/cache/precomputed_quantities.jl:        TD.PhaseNonEquil_pθq(thermo_params, p, θ, q_pt)
src/cache/precomputed_quantities.jl:    return get_ts(ρ, p, θ, e_int, q_tot, q_pt)
src/cache/prognostic_edmf_precomputed_quantities.jl:            TD.dry_pottemp(thermo_params, ᶜts⁰),                               # θ_sat
src/cache/prognostic_edmf_precomputed_quantities.jl:            TD.liquid_ice_pottemp(thermo_params, ᶜts⁰),                        # θ_liq_ice_sat
src/cache/prognostic_edmf_precomputed_quantities.jl:            ),                                                                 # ∂θv∂z_unsat
src/cache/prognostic_edmf_precomputed_quantities.jl:            ),                                                                 # ∂θl∂z_sat
src/initial_conditions/initial_conditions.jl:        θ₀ = FT(280.0)
src/initial_conditions/initial_conditions.jl:        θ = θ₀ * exp(buoy_freq^2 * z / g)
src/initial_conditions/initial_conditions.jl:            g^2 / (cp_d * θ₀ * buoy_freq^2) * (exp(-buoy_freq^2 * z / g) - 1)
src/initial_conditions/initial_conditions.jl:        T = π_exner * θ # temperature
src/initial_conditions/initial_conditions.jl:        θ_b = FT(300)
src/initial_conditions/initial_conditions.jl:        θ_c = FT(-15)
src/initial_conditions/initial_conditions.jl:        θ_p =
src/initial_conditions/initial_conditions.jl:            sqrt(r²) < r_c ? FT(1 / 2) * θ_c * (FT(1) + cospi(sqrt(r²) / r_c)) :
src/initial_conditions/initial_conditions.jl:        θ = θ_b + θ_p # potential temperature
src/initial_conditions/initial_conditions.jl:        π_exn = FT(1) - grav * z / cp_d / θ # exner function
src/initial_conditions/initial_conditions.jl:        T = π_exn * θ # temperature
src/initial_conditions/initial_conditions.jl:        θ_b = FT(300)
src/initial_conditions/initial_conditions.jl:        θ_c = FT(0.5)
src/initial_conditions/initial_conditions.jl:        θ_p =
src/initial_conditions/initial_conditions.jl:            sqrt(r²) < r_c ? FT(1 / 2) * θ_c * (FT(1) + cospi(sqrt(r²) / r_c)) :
src/initial_conditions/initial_conditions.jl:        θ = θ_b + θ_p # potential temperature
src/initial_conditions/initial_conditions.jl:        π_exn = FT(1) - grav * z / cp_d / θ # exner function
src/initial_conditions/initial_conditions.jl:        T = π_exn * θ # temperature
src/initial_conditions/initial_conditions.jl:    hydrostatic_pressure_profile(; thermo_params, p_0, [T, θ, q_tot, z_max])
src/initial_conditions/initial_conditions.jl:given `p(0)`, either `T(z)` or `θ(z)`, and optionally also `q_tot(z)`. If
src/initial_conditions/initial_conditions.jl:the specified profiles T(z), θ(z), and/or q_tot(z) are valid.
src/initial_conditions/initial_conditions.jl:    θ = nothing,
src/initial_conditions/initial_conditions.jl:    ts(p, z, ::Nothing, ::Nothing, _) = error("Either T or θ must be specified")
src/initial_conditions/initial_conditions.jl:    ts(p, z, T::FunctionOrSpline, θ::FunctionOrSpline, _) =
src/initial_conditions/initial_conditions.jl:        error("Only one of T and θ can be specified")
src/initial_conditions/initial_conditions.jl:    ts(p, z, ::Nothing, θ::FunctionOrSpline, ::Nothing) =
src/initial_conditions/initial_conditions.jl:        TD.PhaseDry_pθ(thermo_params, p, oftype(p, θ(z)))
src/initial_conditions/initial_conditions.jl:    ts(p, z, ::Nothing, θ::FunctionOrSpline, q_tot::FunctionOrSpline) =
src/initial_conditions/initial_conditions.jl:        TD.PhaseEquil_pθq(
src/initial_conditions/initial_conditions.jl:            oftype(p, θ(z)),
src/initial_conditions/initial_conditions.jl:    dp_dz(p, z) = -grav * TD.air_density(thermo_params, ts(p, z, T, θ, q_tot))
src/initial_conditions/initial_conditions.jl:    θ_func_name = Symbol(IC, :_θ_liq_ice)
src/initial_conditions/initial_conditions.jl:        θ = APL.$θ_func_name(FT)
src/initial_conditions/initial_conditions.jl:        p = hydrostatic_pressure_profile(; thermo_params, p_0, θ)
src/initial_conditions/initial_conditions.jl:                thermo_state = TD.PhaseDry_pθ(thermo_params, p(z), θ(z)),
src/initial_conditions/initial_conditions.jl:    θ_func_name = Symbol(IC, :_θ_liq_ice)
src/initial_conditions/initial_conditions.jl:        θ = APL.$θ_func_name(FT)
src/initial_conditions/initial_conditions.jl:        p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
src/initial_conditions/initial_conditions.jl:                thermo_state = TD.PhaseEquil_pθq(
src/initial_conditions/initial_conditions.jl:                    θ(z),
src/initial_conditions/initial_conditions.jl:    θ_func_name = Symbol(IC, :_θ_liq_ice)
src/initial_conditions/initial_conditions.jl:        θ = APL.$θ_func_name(FT)
src/initial_conditions/initial_conditions.jl:        p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
src/initial_conditions/initial_conditions.jl:                thermo_state = TD.PhaseEquil_pθq(
src/initial_conditions/initial_conditions.jl:                    θ(z),
src/initial_conditions/initial_conditions.jl:    θ = APL.Rico_θ_liq_ice(FT)
src/initial_conditions/initial_conditions.jl:    p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
src/initial_conditions/initial_conditions.jl:            thermo_state = TD.PhaseEquil_pθq(
src/initial_conditions/initial_conditions.jl:                θ(z),
src/parameterized_tendencies/radiation/held_suarez.jl:    Δθ_z = FT(CAP.Δθ_z(params))
src/parameterized_tendencies/radiation/held_suarez.jl:                    Δθ_z *
src/parameterized_tendencies/radiation/held_suarez.jl:    if :ρθ in propertynames(Y.c)
src/parameterized_tendencies/radiation/held_suarez.jl:        @. Yₜ.c.ρθ[colidx] -=
src/parameterized_tendencies/radiation/radiation.jl:    if :ρθ in propertynames(Y.c)
src/parameterized_tendencies/radiation/radiation.jl:        error("radiation_tendency! not implemented for ρθ")
src/parameters/Parameters.jl:    Δθ_z::FT
src/prognostic_equations/implicit/implicit_tendency.jl:            χ_name == :θ ? energy_upwinding : tracer_upwinding,
src/prognostic_equations/implicit/implicit_solver.jl:        MatrixFields.unrolled_findonly(is_in_Y, (@name(c.ρe_tot), @name(c.ρθ)))
src/prognostic_equations/implicit/implicit_solver.jl:    # ᶜp = p_ref_theta * (ᶜρθ * R_d / p_ref_theta)^(1 / (1 - κ_d)) or
src/prognostic_equations/implicit/implicit_solver.jl:    # ∂(ᶜp)/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #         R_d / (1 - κ_d) * (ᶜρθ * R_d / p_ref_theta)^(κ_d / (1 - κ_d))
src/prognostic_equations/implicit/implicit_solver.jl:        (@name(c.ρθ), @name(ᶜspecific.θ), energy_upwinding),
src/prognostic_equations/implicit/implicit_solver.jl:    # ∂(ᶠu₃ₜ)/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #     ∂(ᶠu₃ₜ)/∂(ᶠgradᵥ(ᶜp - ᶜp_ref)) ⋅ ∂(ᶠgradᵥ(ᶜp - ᶜp_ref))/∂(ᶜρθ) =
src/prognostic_equations/implicit/implicit_solver.jl:    #     ∂(ᶜp)/∂(ᶜρθ) and
src/prognostic_equations/implicit/implicit_solver.jl:    if MatrixFields.has_field(Y, @name(c.ρθ))
src/prognostic_equations/implicit/implicit_solver.jl:        ᶜρθ = Y.c.ρθ
src/prognostic_equations/implicit/implicit_solver.jl:        ∂ᶠu₃_err_∂ᶜρθ = matrix[@name(f.u₃), @name(c.ρθ)]
src/prognostic_equations/implicit/implicit_solver.jl:        @. ∂ᶠu₃_err_∂ᶜρθ[colidx] =
src/prognostic_equations/implicit/implicit_solver.jl:                (ᶜρθ[colidx] * R_d / p_ref_theta)^(κ_d / (1 - κ_d)),
src/prognostic_equations/buoyancy_gradients.jl:    ∂b∂θv = g * (R_d * bg_model.ρ / bg_model.p) * Π
src/prognostic_equations/buoyancy_gradients.jl:            TD.PhaseDry_pθ(thermo_params, bg_model.p, bg_model.θ_liq_ice_sat)
src/prognostic_equations/buoyancy_gradients.jl:            TD.PhaseEquil_pθq(
src/prognostic_equations/buoyancy_gradients.jl:                bg_model.θ_liq_ice_sat,
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂θl_sat = (
src/prognostic_equations/buoyancy_gradients.jl:            ∂b∂θv * (
src/prognostic_equations/buoyancy_gradients.jl:            (lh / cp_m / bg_model.t_sat * ∂b∂θl_sat - ∂b∂θv) * bg_model.θ_sat
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂θl_sat = FT(0)
src/prognostic_equations/buoyancy_gradients.jl:    ∂b∂z = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_sat, ∂b∂qt_sat)
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂θv::FT,
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂θl_sat::FT,
src/prognostic_equations/buoyancy_gradients.jl:    ∂b∂θv::FT,
src/prognostic_equations/buoyancy_gradients.jl:    ∂b∂θl_sat::FT,
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂z_θl_sat = ∂b∂θl_sat * bg_model.∂θl∂z_sat
src/prognostic_equations/buoyancy_gradients.jl:        ∂b∂z_θl_sat = FT(0)
src/prognostic_equations/buoyancy_gradients.jl:        bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_unsat : FT(0)
src/prognostic_equations/buoyancy_gradients.jl:    ∂b∂z_sat = ∂b∂z_θl_sat + ∂b∂z_qt_sat
src/prognostic_equations/cloud_fraction.jl:    θ_liq_ice = TD.liquid_ice_pottemp(thermo_params, ᶜts)
src/prognostic_equations/cloud_fraction.jl:    θl′θl′ = FT(5)
src/prognostic_equations/cloud_fraction.jl:    θl′qt′ = FT(0)
src/prognostic_equations/cloud_fraction.jl:        θ_liq_ice,
src/prognostic_equations/cloud_fraction.jl:        θl′θl′,
src/prognostic_equations/cloud_fraction.jl:        θl′qt′,
src/prognostic_equations/cloud_fraction.jl:    θ_liq_ice,
src/prognostic_equations/cloud_fraction.jl:    θl′θl′,
src/prognostic_equations/cloud_fraction.jl:    θl′qt′,
src/prognostic_equations/cloud_fraction.jl:        θl_mean = θ_liq_ice,
src/prognostic_equations/cloud_fraction.jl:        θl′θl′,
src/prognostic_equations/cloud_fraction.jl:        θl′qt′,
src/prognostic_equations/cloud_fraction.jl:    # θl - liquid ice potential temperature
src/prognostic_equations/cloud_fraction.jl:    (; qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, p_c) = vars
src/prognostic_equations/cloud_fraction.jl:    eps_θ = FT(eps(FT))
src/prognostic_equations/cloud_fraction.jl:    σ_θ::FT = sqrt(θl′θl′)
src/prognostic_equations/cloud_fraction.jl:    _corr::FT = (θl′qt′ / max(σ_q, eps_q))
src/prognostic_equations/cloud_fraction.jl:    corr::FT = max(min(_corr / max(σ_θ, eps_θ), 1), -1)
src/prognostic_equations/cloud_fraction.jl:    σ_c = sqrt(max(1 - corr * corr, 0)) * σ_θ
src/prognostic_equations/cloud_fraction.jl:        μ_c = θl_mean + sqrt2 * corr * σ_θ * χ[1]
src/prognostic_equations/cloud_fraction.jl:        θ_hat = μ_c + sqrt2 * σ_c * χ[2]
src/prognostic_equations/cloud_fraction.jl:        return (θ_hat, q_tot_hat)
src/prognostic_equations/cloud_fraction.jl:        thermo_state(thermo_params; p = p_c, θ = x_hat[1], q_tot = x_hat[2])
src/prognostic_equations/advection.jl:    if :ρθ in propertynames(Y.c)
src/prognostic_equations/advection.jl:        @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu)
src/prognostic_equations/hyperdiffusion.jl:    if :θ in propertynames(ᶜspecific)
src/prognostic_equations/hyperdiffusion.jl:        @. ᶜ∇²specific_energy = wdivₕ(gradₕ(ᶜspecific.θ))
src/prognostic_equations/hyperdiffusion.jl:    ᶜρ_energyₜ = :θ in propertynames(ᶜspecific) ? Yₜ.c.ρθ : Yₜ.c.ρe_tot
src/prognostic_equations/vertical_diffusion_boundary_layer.jl:        elseif χ_name == :θ
src/prognostic_equations/vertical_diffusion_boundary_layer.jl:            @. ρ_flux_χ[colidx] = sfc_conditions.ρ_flux_θ[colidx]
src/solver/types.jl:    θ_sat::FT
src/solver/types.jl:    θ_liq_ice_sat::FT
src/solver/types.jl:    ∂θv∂z_unsat::FT
src/solver/types.jl:    ∂θl∂z_sat::FT
src/surface_conditions/surface_state.jl:    θAndQFluxes(; θ_flux, q_flux)
src/surface_conditions/surface_state.jl:Base.@kwdef struct θAndQFluxes{FT, FTN <: Union{FT, Nothing}} <:
src/surface_conditions/surface_state.jl:    θ_flux::FT
src/surface_conditions/surface_setups.jl:    θ_flux = FT(0.06)
src/surface_conditions/surface_setups.jl:    parameterization = MoninObukhov(; z0, fluxes = θAndQFluxes(; θ_flux))
src/surface_conditions/surface_setups.jl:    θ_flux = FT(0.06)
src/surface_conditions/surface_setups.jl:        MoninObukhov(; z0, fluxes = θAndQFluxes(; θ_flux, q_flux))
src/surface_conditions/surface_setups.jl:    θ_flux = FT(8e-3)
src/surface_conditions/surface_setups.jl:    fluxes = θAndQFluxes(; θ_flux, q_flux)
src/surface_conditions/surface_setups.jl:    θ_flux0 = FT(8e-3)
src/surface_conditions/surface_setups.jl:            θAndQFluxes(; θ_flux = θ_flux0 * weight, q_flux = q_flux0 * weight)
src/surface_conditions/surface_setups.jl:    θ = FT(299)
src/surface_conditions/surface_setups.jl:    ts = TD.PhaseNonEquil_pθq(thermo_params, p, θ, TD.PhasePartition(q_vap))
src/surface_conditions/surface_conditions.jl:            elseif parameterization.fluxes isa θAndQFluxes
src/surface_conditions/surface_conditions.jl:                (; θ_flux, q_flux) = parameterization.fluxes
src/surface_conditions/surface_conditions.jl:                shf = θ_flux * ρ * TD.cp_m(thermo_params, ts)
src/utils/variable_manipulations.jl:Creates a copy of `specific_state` with the energy variable (`θ` or `e_tot`)
src/utils/variable_manipulations.jl:    Base.structdiff(specific_state, NamedTuple{(:θ, :e_tot)})
src/utils/utilities.jl:is_energy_var(symbol) = symbol in (:ρθ, :ρe_tot, :ρaθ, :ρae_tot)

[szy21]

Thanks @Sbozzolo! I'll remove ρθ or θ in the tendencies. Most θ in initial conditions, surface conditions, parameterizations, and plotting scripts should stay though, as they either follow the literature (prescribed θ or θ_flux) or is something we want to have in the plots.

[akshaysridhar]

There are also some remaining references to PotentialTemperature

./post_processing/post_processing_funcs.jl:14:    PotentialTemperature,
./src/solver/types.jl:10:struct PotentialTemperature <: AbstractEnergyFormulation end
./src/surface_conditions/SurfaceConditions.jl:5:import ..PotentialTemperature
./src/initial_conditions/InitialConditions.jl:4:import ..PotentialTemperature

θ_flux or θ doesn't appear to be tied to this abstract type - it seems like we should be able to safely remove references to the PotentialTemperature type while retaining initial condition descriptions based on some θ profile.

[Sbozzolo]

Thanks @Sbozzolo! I'll remove ρθ or θ in the tendencies. Most θ in initial conditions, surface conditions, parameterizations, and plotting scripts should stay though, as they either follow the literature (prescribed θ or θ_flux) or is something we want to have in the plots.

Fina, as long as they don't depend explicitly on the removed variable

[szy21]

@Sbozzolo I think I have removed everything that depends on the removed variable. Would you like to take a look again?

[Sbozzolo]

@Sbozzolo I think I have removed everything that depends on the removed variable. Would you like to take a look again?

Yes, I'll do it now.

[akshaysridhar]

There are also some remaining references to PotentialTemperature

./post_processing/post_processing_funcs.jl:14:    PotentialTemperature,
./src/solver/types.jl:10:struct PotentialTemperature <: AbstractEnergyFormulation end
./src/surface_conditions/SurfaceConditions.jl:5:import ..PotentialTemperature
./src/initial_conditions/InitialConditions.jl:4:import ..PotentialTemperature

θ_flux or θ doesn't appear to be tied to this abstract type - it seems like we should be able to safely remove references to the PotentialTemperature type while retaining initial condition descriptions based on some θ profile.

This is now resolved as of commit af6fad0 . Thanks for making these changes @szy21 !

[Sbozzolo]

The function thermo_state in precomputed_quantities.jl can probably be simplified given that theta is no longer used.

I think it'd look like:

function thermo_state(
    thermo_params;
    ρ = nothing,
    p = nothing,
    e_int = nothing,
    q_tot = nothing,
    q_pt = nothing,
)
    get_ts(ρ::Real, ::Nothing, e_int::Real, ::Nothing, ::Nothing) =
        TD.PhaseDry_ρe(thermo_params, ρ, e_int)
    get_ts(ρ::Real, ::Nothing, e_int::Real, q_tot::Real, ::Nothing) =
        TD.PhaseEquil_ρeq(
            thermo_params,
            ρ,
            e_int,
            q_tot,
            3,
            eltype(thermo_params)(0.003),
        )
    get_ts(ρ::Real, ::Nothing, e_int::Real, ::Nothing, q_pt) =
        TD.PhaseNonEquil_ρeq(thermo_params, ρ, e_int, q_pt)
    get_ts(::Nothing, p::Real, e_int::Real, ::Nothing, ::Nothing) =
        TD.PhaseDry_pe(thermo_params, p, e_int)
    get_ts(::Nothing, p::Real, e_int::Real, q_tot::Real, ::Nothing) =
        TD.PhaseEquil_peq(thermo_params, p, e_int, q_tot)
    get_ts(::Nothing, p::Real, e_int::Real, ::Nothing, q_pt) =
        TD.PhaseNonEquil_peq(thermo_params, p, e_int, q_pt)
    return get_ts(ρ, p, e_int, q_tot, q_pt)
end

[szy21]

We still use θ to construct thermo state sometimes, for example in initial conditions here and in closures here. This is because some literature prescribe θ for initial conditions, and some of our closures is dependent on θ (as we haven't developed the same closure for e_tot yet).

[Sbozzolo]

We still use θ to construct thermo state sometimes, for example in initial conditions here and in closures here. This is because some literature prescribe θ for initial conditions, and some of our closures is dependent on θ (as we haven't developed the same closure for e_tot yet).

Okay, thank you!

Everything else is looking good.