diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index fe8e81d4..6c37bb6c 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -62,5 +62,5 @@ steps:
       slurm_gres: "gpu:1"
     env:
       JULIA_VERSION: "1.7.0"
-      CUDA_VERSION: "11.2"
+      CUDA_VERSION: "10.2"
     timeout_in_minutes: 60
diff --git a/Project.toml b/Project.toml
index 291532f4..460a84fc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,12 +1,11 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.5.8"
+version = "0.5.9"
 
 [deps]
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
@@ -16,5 +15,4 @@ CLIMAParameters = "0.1, 0.2, 0.3"
 DocStringExtensions = "0.8.1"
 KernelAbstractions = "0.5, 0.6, 0.7"
 RootSolvers = "0.2"
-ExprTools = "0.1 - 0.1.3"
 julia = "1.4"
diff --git a/perf/alloc_per_constructor.jl b/perf/alloc_per_constructor.jl
index 44d8700e..af98187c 100644
--- a/perf/alloc_per_constructor.jl
+++ b/perf/alloc_per_constructor.jl
@@ -8,15 +8,12 @@ include("common.jl")
 constructor = ENV["ALLOCATION_CONSTRUCTOR"]
 @info "Recording allocations for $constructor"
 
-
-if constructor == "ρeq"
-    thermo_state = thermo_state_ρeq_newton
-else
-end
-if constructor == "pθq"
-    thermo_state = thermo_state_pθq
-else
-end
+thermo_state_map = Dict(
+    "ρeq" => thermo_state_ρeq_newton,
+    "pθq" => thermo_state_pθq,
+    "pTq" => thermo_state_pTq,
+)
+thermo_state = thermo_state_map[constructor]
 
 thermo_state() # compile first
 Profile.clear_malloc_data()
diff --git a/perf/allocs.jl b/perf/allocs.jl
index 85d35c55..84264b8c 100644
--- a/perf/allocs.jl
+++ b/perf/allocs.jl
@@ -20,7 +20,7 @@ all_dirs_to_monitor = [
 # https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
 ENV["GKSwstype"] = "nul"
 
-constructors = ["ρeq", "pθq"]
+constructors = ["ρeq", "pθq", "pTq"]
 
 allocs = Dict()
 for constructor in constructors
@@ -101,6 +101,10 @@ function plot_allocs(constructor, allocs_per_case, n_unique_bytes)
     end
 
     all_bytes = all_bytes ./ 10^3
+    if isempty(all_bytes)
+        @info "$constructor: 0 allocations! 🎉"
+        return nothing
+    end
     max_bytes = maximum(all_bytes)
     @info "$constructor: $all_bytes"
     xtick_name(filename, linenumber) = "$filename, line number: $linenumber"
diff --git a/perf/benchmark.jl b/perf/benchmark.jl
index 3ebb3fb9..cb7877bd 100644
--- a/perf/benchmark.jl
+++ b/perf/benchmark.jl
@@ -10,5 +10,7 @@ include("common.jl")
     show(stdout, MIME("text/plain"), trial)
     trial = BenchmarkTools.@benchmark $thermo_state_pθq()
     show(stdout, MIME("text/plain"), trial)
+    trial = BenchmarkTools.@benchmark $thermo_state_pTq()
+    show(stdout, MIME("text/plain"), trial)
 
 end
diff --git a/perf/common.jl b/perf/common.jl
index 841ea8b3..8ccf9ae3 100644
--- a/perf/common.jl
+++ b/perf/common.jl
@@ -19,7 +19,7 @@ const param_set = EarthParameterSet()
 
 ArrayType = Array{Float64}
 profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
-UnPack.@unpack e_int, ρ, p, θ_liq_ice, q_tot = profiles
+UnPack.@unpack e_int, T, ρ, p, θ_liq_ice, q_tot = profiles
 
 function thermo_state_ρeq_newton()
     ts =
@@ -57,3 +57,7 @@ end
 function thermo_state_pθq()
     ts = TD.PhaseEquil_pθq.(param_set, p, θ_liq_ice, q_tot)
 end
+
+function thermo_state_pTq()
+    ts = TD.PhaseEquil_pTq.(param_set, p, T, q_tot)
+end
diff --git a/src/TestedProfiles.jl b/src/TestedProfiles.jl
index 8ebd54e8..0283e5b9 100644
--- a/src/TestedProfiles.jl
+++ b/src/TestedProfiles.jl
@@ -161,7 +161,7 @@ end
 Returns a `ProfileSet` used to test dry thermodynamic states.
 """
 function PhaseDryProfiles(param_set::APS, ::Type{ArrayType}) where {ArrayType}
-    phase_type = TD.PhaseDry
+    phase_type = TD.PhaseDry{eltype(ArrayType), typeof(param_set)}
 
     z_range, relative_sat, T_surface, T_min = input_config(ArrayType)
     z, T_virt, p, RS =
@@ -221,7 +221,7 @@ end
 Returns a `ProfileSet` used to test moist states in thermodynamic equilibrium.
 """
 function PhaseEquilProfiles(param_set::APS, ::Type{ArrayType}) where {ArrayType}
-    phase_type = TD.PhaseEquil
+    phase_type = TD.PhaseEquil{eltype(ArrayType), typeof(param_set)}
 
     # Prescribe z_range, relative_sat, T_surface, T_min
     z_range, relative_sat, T_surface, T_min = input_config(ArrayType)
diff --git a/src/config_numerical_method.jl b/src/config_numerical_method.jl
index cbcbe5ee..868bf7ac 100644
--- a/src/config_numerical_method.jl
+++ b/src/config_numerical_method.jl
@@ -30,8 +30,8 @@ function sa_numerical_method(
     ρ::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.NewtonsMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.NewtonsMethod, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     T_init =
         max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
@@ -47,8 +47,8 @@ function sa_numerical_method(
     ρ::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.NewtonsMethodAD}
+    ::Type{phase_type},
+) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     T_init =
         max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
@@ -61,8 +61,8 @@ function sa_numerical_method(
     ρ::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.SecantMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
     T_2 = air_temperature(param_set, e_int, q_pt)
@@ -77,8 +77,8 @@ function sa_numerical_method(
     ρ::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.RegulaFalsiMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
     T_2 = air_temperature(param_set, e_int, q_pt)
@@ -97,8 +97,8 @@ function sa_numerical_method_ρpq(
     ρ::FT,
     p::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.NewtonsMethodAD}
+    ::Type{phase_type},
+) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil}
     q_pt = PhasePartition(q_tot)
     T_init = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt)
     return RS.NewtonsMethodAD(T_init)
@@ -110,8 +110,8 @@ function sa_numerical_method_ρpq(
     ρ::FT,
     p::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.RegulaFalsiMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil}
     q_pt = PhasePartition(q_tot)
     T_1 = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) - 5
     T_2 = air_temperature_from_ideal_gas_law(param_set, p, ρ, q_pt) + 5
@@ -128,8 +128,8 @@ function sa_numerical_method_peq(
     p::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.NewtonsMethodAD}
+    ::Type{phase_type},
+) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     T_init =
         max(T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor
@@ -142,8 +142,8 @@ function sa_numerical_method_peq(
     p::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.SecantMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     q_pt = PhasePartition(q_tot, FT(0), q_tot) # Assume all ice
     T_2 = air_temperature(param_set, e_int, q_pt)
@@ -162,8 +162,8 @@ function sa_numerical_method_pθq(
     p::FT,
     θ_liq_ice::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.RegulaFalsiMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.RegulaFalsiMethod, phase_type <: PhaseEquil}
     _T_min::FT = CPP.T_min(param_set)
     _T_max::FT = CPP.T_max(param_set)
     air_temp(q) = air_temperature_given_pθq(param_set, p, θ_liq_ice, q)
@@ -179,8 +179,8 @@ function sa_numerical_method_pθq(
     p::FT,
     θ_liq_ice::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.SecantMethod}
+    ::Type{phase_type},
+) where {FT, NM <: RS.SecantMethod, phase_type <: PhaseEquil}
     _T_min::FT = CPP.T_min(param_set)
     air_temp(q) = air_temperature_given_pθq(param_set, p, θ_liq_ice, q)
     T_1 = max(_T_min, air_temp(PhasePartition(q_tot))) # Assume all vapor
@@ -195,8 +195,8 @@ function sa_numerical_method_pθq(
     p::FT,
     θ_liq_ice::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT, NM <: RS.NewtonsMethodAD}
+    ::Type{phase_type},
+) where {FT, NM <: RS.NewtonsMethodAD, phase_type <: PhaseEquil}
     T_min::FT = CPP.T_min(param_set)
     air_temp(q) = air_temperature_given_pθq(param_set, p, θ_liq_ice, q)
     T_init = max(T_min, air_temp(PhasePartition(q_tot))) # Assume all vapor
diff --git a/src/relations.jl b/src/relations.jl
index ccea2673..5d21dc88 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -526,8 +526,8 @@ function internal_energy_sat(
     T::FT,
     ρ::FT,
     q_tot::FT,
-    phase_type::Type{<:ThermodynamicState},
-) where {FT <: Real}
+    ::Type{phase_type},
+) where {FT <: Real, phase_type <: ThermodynamicState}
     return internal_energy(
         param_set,
         T,
@@ -851,9 +851,9 @@ end
 
 function saturation_vapor_pressure(
     param_set::APS,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     T::FT,
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     LH_s0 = FT(CPP.LH_s0(param_set))
     LH_v0 = FT(CPP.LH_v0(param_set))
     cp_v = FT(CPP.cp_v(param_set))
@@ -933,9 +933,9 @@ function q_vap_saturation(
     param_set::APS,
     T::FT,
     ρ::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q::PhasePartition{FT} = q_pt_0(FT),
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
 
     LH_v0::FT = CPP.LH_v0(param_set)
     LH_s0::FT = CPP.LH_s0(param_set)
@@ -1032,8 +1032,8 @@ function q_vap_saturation_from_pressure(
     q_tot::FT,
     p::FT,
     T::FT,
-    phase_type::Type{<:ThermodynamicState},
-) where {FT <: Real}
+    ::Type{phase_type},
+) where {FT <: Real, phase_type <: ThermodynamicState}
     R_v::FT = CPP.R_v(param_set)
     R_d::FT = CPP.R_d(param_set)
     p_v_sat = saturation_vapor_pressure(param_set, phase_type, T)
@@ -1108,9 +1108,9 @@ function saturation_excess(
     param_set::APS,
     T::FT,
     ρ::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q::PhasePartition{FT},
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     return max(0, q.tot - q_vap_saturation(param_set, T, ρ, phase_type, q))
 end
 
@@ -1168,9 +1168,9 @@ is a function that is 1 above `T_freeze` and goes to zero below `T_freeze`.
 function liquid_fraction(
     param_set::APS,
     T::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q::PhasePartition{FT} = q_pt_0(FT),
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     _T_freeze::FT = CPP.T_freeze(param_set)
     return FT(T > _T_freeze)
 end
@@ -1178,9 +1178,9 @@ end
 function liquid_fraction(
     param_set::APS,
     T::FT,
-    phase_type::Type{<:PhaseNonEquil},
+    ::Type{phase_type},
     q::PhasePartition{FT} = q_pt_0(FT),
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: PhaseNonEquil}
     q_c = condensate(q)     # condensate specific humidity
     if has_condensate(q_c)
         return q.liq / q_c
@@ -1220,8 +1220,8 @@ function PhasePartition_equil(
     T::FT,
     ρ::FT,
     q_tot::FT,
-    phase_type::Type{<:ThermodynamicState},
-) where {FT <: Real}
+    ::Type{phase_type},
+) where {FT <: Real, phase_type <: ThermodynamicState}
     _liquid_frac = liquid_fraction(param_set, T, phase_type)                    # fraction of condensate that is liquid
     q_c = saturation_excess(param_set, T, ρ, phase_type, PhasePartition(q_tot)) # condensate specific humidity
     q_liq = _liquid_frac * q_c                                                  # liquid specific humidity
@@ -1255,8 +1255,8 @@ function PhasePartition_equil_given_p(
     T::FT,
     p::FT,
     q_tot::FT,
-    phase_type::Type{<:ThermodynamicState},
-) where {FT <: Real}
+    ::Type{phase_type},
+) where {FT <: Real, phase_type <: ThermodynamicState}
 
     q_v_sat = q_vap_saturation_from_pressure(param_set, q_tot, p, T, phase_type)
     _liquid_frac = liquid_fraction(param_set, T, phase_type)
@@ -1282,8 +1282,8 @@ function ∂e_int_∂T(
     e_int::FT,
     ρ::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
-) where {FT <: Real}
+    ::Type{phase_type},
+) where {FT <: Real, phase_type <: PhaseEquil}
     LH_v0::FT = CPP.LH_v0(param_set)
     LH_s0::FT = CPP.LH_s0(param_set)
     R_v::FT = CPP.R_v(param_set)
@@ -1346,10 +1346,10 @@ function saturation_adjustment(
     e_int::FT,
     ρ::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
+    ::Type{phase_type},
     maxiter::Int,
     temperature_tol::FT,
-) where {FT <: Real, sat_adjust_method}
+) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil}
     _T_min::FT = CPP.T_min(param_set)
     cv_d::FT = CPP.cv_d(param_set)
     # Convert temperature tolerance to a convergence criterion on internal energy residuals
@@ -1442,10 +1442,10 @@ function saturation_adjustment_given_peq(
     p::FT,
     e_int::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
+    ::Type{phase_type},
     maxiter::Int,
     temperature_tol::FT,
-) where {FT <: Real, sat_adjust_method}
+) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil}
     _T_min::FT = CPP.T_min(param_set)
     cv_d = FT(CPP.cv_d(param_set))
     # Convert temperature tolerance to a convergence criterion on internal energy residuals
@@ -1541,10 +1541,10 @@ function saturation_adjustment_ρpq(
     ρ::FT,
     p::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
+    ::Type{phase_type},
     maxiter::Int,
     temperature_tol::FT = sqrt(eps(FT)),
-) where {FT <: Real, sat_adjust_method}
+) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil}
     tol = RS.SolutionTolerance(temperature_tol)
     # Use `oftype` to preserve diagonalized type signatures:
     sol = RS.find_zero(
@@ -1651,10 +1651,10 @@ function saturation_adjustment_given_ρθq(
     ρ::FT,
     θ_liq_ice::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
+    ::Type{phase_type},
     maxiter::Int,
     tol::RS.AbstractTolerance,
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: PhaseEquil}
     _T_min::FT = CPP.T_min(param_set)
     air_temp(q) = air_temperature_given_ρθq(param_set, ρ, θ_liq_ice, q)
     T_1 = max(_T_min, air_temp(PhasePartition(q_tot))) # Assume all vapor
@@ -1733,10 +1733,10 @@ function saturation_adjustment_given_pθq(
     p::FT,
     θ_liq_ice::FT,
     q_tot::FT,
-    phase_type::Type{<:PhaseEquil},
+    ::Type{phase_type},
     maxiter::Int,
     temperature_tol::FT,
-) where {FT <: Real, sat_adjust_method}
+) where {FT <: Real, sat_adjust_method, phase_type <: PhaseEquil}
     tol = RS.ResidualTolerance(temperature_tol)
     T_min::FT = CPP.T_min(param_set)
     T_freeze::FT = CPP.T_freeze(param_set)
@@ -1940,8 +1940,8 @@ function virt_temp_from_RH(
     T::FT,
     ρ::FT,
     RH::FT,
-    phase_type::Type{<:ThermodynamicState},
-) where {FT <: AbstractFloat}
+    ::Type{phase_type},
+) where {FT <: AbstractFloat, phase_type <: ThermodynamicState}
     q_tot = RH * q_vap_saturation(param_set, T, ρ, phase_type)
     q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
     return virtual_temperature(param_set, T, ρ, q_pt)
@@ -1962,10 +1962,10 @@ function temperature_and_humidity_given_TᵥρRH(
     T_virt::FT,
     ρ::FT,
     RH::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     maxiter::Int = 100,
     tol::RS.AbstractTolerance = RS.ResidualTolerance{FT}(sqrt(eps(FT))),
-) where {FT <: AbstractFloat}
+) where {FT <: AbstractFloat, phase_type <: ThermodynamicState}
 
     _T_min::FT = CPP.T_min(param_set)
     _T_max = T_virt
@@ -2210,9 +2210,9 @@ function liquid_ice_pottemp_sat(
     param_set::APS,
     T::FT,
     ρ::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q::PhasePartition{FT} = q_pt_0(FT),
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     q_v_sat = q_vap_saturation(param_set, T, ρ, phase_type, q)
     return liquid_ice_pottemp(param_set, T, ρ, PhasePartition(q_v_sat))
 end
@@ -2232,9 +2232,9 @@ function liquid_ice_pottemp_sat(
     param_set::APS,
     T::FT,
     ρ::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q_tot::FT,
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     return liquid_ice_pottemp(
         param_set,
         T,
@@ -2377,9 +2377,9 @@ function relative_humidity(
     param_set::APS,
     T::FT,
     p::FT,
-    phase_type::Type{<:ThermodynamicState},
+    ::Type{phase_type},
     q::PhasePartition{FT} = q_pt_0(FT),
-) where {FT <: Real}
+) where {FT <: Real, phase_type <: ThermodynamicState}
     R_v::FT = CPP.R_v(param_set)
     q_vap = vapor_specific_humidity(q)
     p_vap = q_vap * air_density(param_set, T, p, q) * R_v * T
diff --git a/src/states.jl b/src/states.jl
index 6b9b302a..fe81d79c 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -253,7 +253,7 @@ function PhaseEquil_ρeq(
 ) where {FT <: Real, sat_adjust_method, IT <: ITERTYPE, FTT <: TOLTYPE(FT)}
     maxiter === nothing && (maxiter = 8)
     temperature_tol === nothing && (temperature_tol = FT(1e-1))
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment(
         sat_adjust_method,
@@ -323,7 +323,7 @@ function PhaseEquil_ρθq(
 ) where {FT <: Real, IT <: ITERTYPE, FTT <: TOLTYPE(FT)}
     maxiter === nothing && (maxiter = 36)
     temperature_tol === nothing && (temperature_tol = FT(1e-1))
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     tol = RS.ResidualTolerance(temperature_tol)
     T = saturation_adjustment_given_ρθq(
         param_set,
@@ -356,7 +356,7 @@ function PhaseEquil_ρTq(
     T::FT,
     q_tot::FT,
 ) where {FT <: Real}
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     q_pt = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
     p = air_pressure(param_set, T, ρ, q_pt)
     e_int = internal_energy(param_set, T, q_pt)
@@ -379,7 +379,7 @@ function PhaseEquil_pTq(
     T::FT,
     q_tot::FT,
 ) where {FT <: Real}
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     q_pt = PhasePartition_equil_given_p(param_set, T, p, q_tot_safe, phase_type)
     ρ = air_density(param_set, T, p, q_pt)
@@ -415,7 +415,7 @@ function PhaseEquil_peq(
 ) where {FT <: Real, sat_adjust_method, IT <: ITERTYPE, FTT <: TOLTYPE(FT)}
     maxiter === nothing && (maxiter = 40)
     temperature_tol === nothing && (temperature_tol = FT(1e-2))
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment_given_peq(
         sat_adjust_method,
@@ -463,7 +463,7 @@ function PhaseEquil_ρpq(
     ::Type{sat_adjust_method} = RS.NewtonsMethodAD,
 ) where {FT <: Real, sat_adjust_method}
 
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     if perform_sat_adjust
         T = saturation_adjustment_ρpq(
             sat_adjust_method,
@@ -509,7 +509,7 @@ function PhaseEquil_pθq(
 ) where {FT <: Real, IT <: ITERTYPE, FTT <: TOLTYPE(FT), sat_adjust_method}
     maxiter === nothing && (maxiter = 50)
     temperature_tol === nothing && (temperature_tol = FT(1e-3))
-    phase_type = PhaseEquil
+    phase_type = PhaseEquil{FT, typeof(param_set)}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment_given_pθq(
         sat_adjust_method,
@@ -614,7 +614,7 @@ function PhaseNonEquil_ρθq(
     maxiter::Int = 10,
     potential_temperature_tol::FT = FT(1e-2),
 ) where {FT <: Real}
-    phase_type = PhaseNonEquil
+    phase_type = PhaseNonEquil{FT, typeof(param_set)}
     tol = RS.ResidualTolerance(potential_temperature_tol)
     T = air_temperature_given_ρθq_nonlinear(
         param_set,