diff --git a/docs/src/saturation_adjustment.jl b/docs/src/saturation_adjustment.jl
index 82f3af9f..87948408 100644
--- a/docs/src/saturation_adjustment.jl
+++ b/docs/src/saturation_adjustment.jl
@@ -6,7 +6,6 @@ import Thermodynamics.Parameters as TP
 import CLIMAParameters as CP
 
 FT = Float64
-
 const param_set = TP.ThermodynamicsParameters(FT)
 
 profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
@@ -55,8 +54,8 @@ ts_no_err = Dict(
             )
 
             @inbounds for NM in numerical_methods
-                TD.error_on_non_convergence() = false
                 ts_no_err[NM][n] = TD.PhaseEquil_dev_only(
+                    TD.NullLogger(),
                     param_set,
                     ρ[i],
                     e_int,
@@ -64,10 +63,10 @@ ts_no_err = Dict(
                     sat_adjust_method = NM,
                     maxiter = 10,
                 )
-                TD.error_on_non_convergence() = true
                 # @show n / prod(dims) * 100
                 try
                     ts[NM][n] = TD.PhaseEquil_dev_only(
+                        TD.ErrorLogger(),
                         param_set,
                         ρ[i],
                         e_int,
diff --git a/perf/jet.jl b/perf/jet.jl
index 739a8024..232beede 100644
--- a/perf/jet.jl
+++ b/perf/jet.jl
@@ -1,9 +1,5 @@
 include("common_micro_bm.jl")
 
-# JET does not like KernelAbstractions.@print calls,
-# so we disable them here.
-TD.print_warning() = false
-
 function jet_thermo_states(::Type{FT}) where {FT}
     param_set = TP.ThermodynamicsParameters(FT)
     ArrayType = Array{FT}
@@ -22,7 +18,7 @@ function jet_thermo_states(::Type{FT}) where {FT}
         )
             for cond in conditions(C)
                 args = sample_args(profiles, param_set, cond, C)
-                JET.@test_opt C(param_set, args...)
+                JET.@test_opt C(TD.NullLogger(), param_set, args...)
             end
         end
     end
diff --git a/src/DataCollection.jl b/src/DataCollection.jl
index ae87d46e..a335d515 100644
--- a/src/DataCollection.jl
+++ b/src/DataCollection.jl
@@ -12,10 +12,9 @@ import Thermodynamics as TD
 import RootSolvers as RS
 
 function do_work()
-    # Calls TD.PhaseEquil_ρeq()..., possibly many times
+    # Calls TD.PhaseEquil_ρeq(VerboseLogger(WarningLogger()), )..., possibly many times
 end
 
-TD.solution_type() = RS.VerboseSolution()
 do_work()
 TD.DataCollection.print_summary()
 ```
diff --git a/src/Thermodynamics.jl b/src/Thermodynamics.jl
index 13edeef3..749b8cd2 100644
--- a/src/Thermodynamics.jl
+++ b/src/Thermodynamics.jl
@@ -56,26 +56,14 @@ import .Parameters
 const TP = Parameters
 const APS = TP.ThermodynamicsParameters
 
-# Allow users to skip error on non-convergence
-# by importing:
-# ```julia
-# import Thermodynamics
-# Thermodynamics.error_on_non_convergence() = false
-# ```
-# Error on convergence must be the default
-# behavior because this can result in printing
-# very large logs resulting in CI to seemingly hang.
-@inline error_on_non_convergence() = true
-
-# Allow users to skip printing warnings on non-convergence
-@inline print_warning() = true
-
 @inline q_pt_0(::Type{FT}) where {FT} = PhasePartition(FT(0), FT(0), FT(0))
 
-@inline solution_type() = RS.CompactSolution()
 include("DataCollection.jl")
 import .DataCollection
 
+include("logger.jl")
+solution_type(::AbstractThermodynamicsLogger) = RS.CompactSolution()
+solution_type(::VerboseLogger) = RS.CompactSolution()
 include("states.jl")
 include("relations.jl")
 include("isentropic.jl")
diff --git a/src/logger.jl b/src/logger.jl
new file mode 100644
index 00000000..27283a07
--- /dev/null
+++ b/src/logger.jl
@@ -0,0 +1,63 @@
+"""
+    AbstractThermodynamicsLogger
+
+A Thermodynamics logger can be passed to
+thermo state constructors that (may) perform
+saturation adjustment. This logger will statically
+decide if a warning statement is printed, or
+if the function should error on non-convergence.
+"""
+abstract type AbstractThermodynamicsLogger end
+const ATL = AbstractThermodynamicsLogger
+
+"""
+    WarningLogger
+
+Warn when saturation adjustment did not converge.
+"""
+struct WarningLogger <: AbstractThermodynamicsLogger end
+
+"""
+    ErrorLogger
+
+Error (but do not warn) when saturation adjustment
+did not converge.
+
+This is really only useful in our development docs.
+"""
+struct ErrorLogger <: AbstractThermodynamicsLogger end
+
+"""
+    WarnAndErrorLogger
+
+Warn and error (default) when saturation adjustment did not converge.
+"""
+struct WarnAndErrorLogger <: AbstractThermodynamicsLogger end
+
+"""
+    VerboseLogger
+
+A verbose logger, for when we use RootSolvers `VerboseSolution` type.
+"""
+struct VerboseLogger{L} <: AbstractThermodynamicsLogger
+    logger::L
+end
+
+"""
+    NullLogger
+
+Do nothing when saturation adjustment did not converge.
+"""
+struct NullLogger <: AbstractThermodynamicsLogger end
+
+warn_msg(::WarningLogger) = true
+warn_msg(::WarnAndErrorLogger) = true
+warn_msg(l::VerboseLogger) = warn_msg(l.logger)
+warn_msg(::AbstractThermodynamicsLogger) = false
+
+error_msg(::ErrorLogger) = true
+error_msg(::WarnAndErrorLogger) = true
+error_msg(l::VerboseLogger) = error_msg(l.logger)
+error_msg(::AbstractThermodynamicsLogger) = false
+
+Base.broadcastable(x::AbstractThermodynamicsLogger) = tuple(x)
diff --git a/src/relations.jl b/src/relations.jl
index 051fa174..8ce23dc4 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -1521,6 +1521,7 @@ end
 
 """
     saturation_adjustment(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         e_int,
@@ -1554,6 +1555,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     e_int::FT,
@@ -1631,34 +1633,33 @@ See also [`saturation_adjustment`](@ref).
             phase_type,
             T_guess,
         ),
-        solution_type(),
+        solution_type(logger),
         tol,
         maxiter,
     )
 
     DataCollection.log_meta(sol)
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment:\n")
-            print_numerical_method(sat_adjust_method)
-            print_T_guess(sat_adjust_method, T_guess)
-            KA.@print(", e_int=", e_int)
-            KA.@print(", ρ=", ρ)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment:\n")
+        print_numerical_method(sat_adjust_method)
+        print_T_guess(sat_adjust_method, T_guess)
+        KA.@print(", e_int=", e_int)
+        KA.@print(", ρ=", ρ)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
 
 """
     saturation_adjustment_given_peq(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         e_int,
@@ -1694,6 +1695,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment_given_peq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -1737,26 +1739,24 @@ See also [`saturation_adjustment`](@ref).
             phase_type,
             T_guess,
         ),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment_peq:\n")
-            print_numerical_method(sat_adjust_method)
-            print_T_guess(sat_adjust_method, T_guess)
-            KA.@print(", e_int=", e_int)
-            KA.@print(", p=", p)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment_peq:\n")
+        print_numerical_method(sat_adjust_method)
+        print_T_guess(sat_adjust_method, T_guess)
+        KA.@print(", e_int=", e_int)
+        KA.@print(", p=", p)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
@@ -1764,6 +1764,7 @@ end
 
 """
     saturation_adjustment_given_phq(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         h,
@@ -1799,6 +1800,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment_given_phq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -1850,32 +1852,31 @@ See also [`saturation_adjustment`](@ref).
             phase_type,
             T_guess,
         ),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment_phq:\n")
-            print_numerical_method(sat_adjust_method)
-            print_T_guess(sat_adjust_method, T_guess)
-            KA.@print(", h=", h)
-            KA.@print(", p=", p)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment_phq:\n")
+        print_numerical_method(sat_adjust_method)
+        print_T_guess(sat_adjust_method, T_guess)
+        KA.@print(", h=", h)
+        KA.@print(", p=", p)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
 
 """
     saturation_adjustment_ρpq(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         ρ,
@@ -1911,6 +1912,7 @@ using Newtons method using ForwardDiff.
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment_ρpq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     ρ::FT,
@@ -1947,26 +1949,24 @@ See also [`saturation_adjustment`](@ref).
             phase_type,
             T_guess,
         ),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment_ρpq:\n")
-            print_numerical_method(sat_adjust_method)
-            print_T_guess(sat_adjust_method, T_guess)
-            KA.@print(", ρ=", ρ)
-            KA.@print(", p=", p)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment_ρpq:\n")
+        print_numerical_method(sat_adjust_method)
+        print_T_guess(sat_adjust_method, T_guess)
+        KA.@print(", ρ=", ρ)
+        KA.@print(", p=", p)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
@@ -1998,6 +1998,7 @@ end
 
 """
     saturation_adjustment_given_ρθq(
+        logger::ATL,
         param_set,
         ρ,
         θ_liq_ice,
@@ -2029,6 +2030,7 @@ by finding the root of
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment_given_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -2054,31 +2056,30 @@ See also [`saturation_adjustment`](@ref).
     sol = RS.find_zero(
         roots,
         RS.SecantMethod(T_1, T_2),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment_given_ρθq:\n")
-            KA.@print("    Method=SecantMethod")
-            KA.@print(", ρ=", ρ)
-            KA.@print(", θ_liq_ice=", θ_liq_ice)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment_given_ρθq:\n")
+        KA.@print("    Method=SecantMethod")
+        KA.@print(", ρ=", ρ)
+        KA.@print(", θ_liq_ice=", θ_liq_ice)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
 
 """
     saturation_adjustment_given_pθq(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         p,
@@ -2109,6 +2110,7 @@ by finding the root of
 See also [`saturation_adjustment`](@ref).
 """
 @inline function saturation_adjustment_given_pθq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -2169,26 +2171,24 @@ See also [`saturation_adjustment`](@ref).
             phase_type,
             T_guess,
         ),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print("maxiter reached in saturation_adjustment_given_pθq:\n")
-            print_numerical_method(sat_adjust_method)
-            print_T_guess(sat_adjust_method, T_guess)
-            KA.@print(", p=", p)
-            KA.@print(", θ_liq_ice=", θ_liq_ice)
-            KA.@print(", q_tot=", q_tot)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in saturation_adjustment_given_pθq:\n")
+        print_numerical_method(sat_adjust_method)
+        print_T_guess(sat_adjust_method, T_guess)
+        KA.@print(", p=", p)
+        KA.@print(", θ_liq_ice=", θ_liq_ice)
+        KA.@print(", q_tot=", q_tot)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
@@ -2343,7 +2343,7 @@ The dry potential temperature, given a thermodynamic state `ts`.
     return virtual_temperature(param_set, T, ρ, q_pt)
 end
 """
-    temperature_and_humidity_given_TᵥρRH(param_set, T_virt, ρ, RH)
+    temperature_and_humidity_given_TᵥρRH([logger, ]param_set, T_virt, ρ, RH)
 
 The air temperature and `q_tot` where
 
@@ -2354,6 +2354,7 @@ The air temperature and `q_tot` where
  - `phase_type` a thermodynamic state type
 """
 @inline function temperature_and_humidity_given_TᵥρRH(
+    logger::ATL,
     param_set::APS,
     T_virt::FT,
     ρ::FT,
@@ -2370,27 +2371,23 @@ The air temperature and `q_tot` where
     sol = RS.find_zero(
         roots,
         RS.SecantMethod(_T_min, _T_max),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print(
-                "maxiter reached in temperature_and_humidity_given_TᵥρRH:\n"
-            )
-            KA.@print("    Method=SecantMethod")
-            KA.@print(", T_virt=", T_virt)
-            KA.@print(", RH=", RH)
-            KA.@print(", ρ=", ρ)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in temperature_and_humidity_given_TᵥρRH:\n")
+        KA.@print("    Method=SecantMethod")
+        KA.@print(", T_virt=", T_virt)
+        KA.@print(", RH=", RH)
+        KA.@print(", ρ=", ρ)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     T = sol.root
 
@@ -2431,7 +2428,7 @@ and, optionally,
 end
 
 """
-    air_temperature_given_ρθq_nonlinear(param_set, ρ, θ_liq_ice, q::PhasePartition)
+    air_temperature_given_ρθq_nonlinear([logger, ]param_set, ρ, θ_liq_ice, q::PhasePartition)
 
 Computes temperature `T` given
 
@@ -2452,6 +2449,7 @@ by finding the root of
                                q) = 0`
 """
 @inline function air_temperature_given_ρθq_nonlinear(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -2471,29 +2469,25 @@ by finding the root of
     sol = RS.find_zero(
         roots,
         RS.SecantMethod(_T_min, _T_max),
-        RS.CompactSolution(),
+        solution_type(logger),
         tol,
         maxiter,
     )
-    if !sol.converged
-        if print_warning()
-            KA.@print("-----------------------------------------\n")
-            KA.@print(
-                "maxiter reached in air_temperature_given_ρθq_nonlinear:\n"
-            )
-            KA.@print("    Method=SecantMethod")
-            KA.@print(", θ_liq_ice=", θ_liq_ice)
-            KA.@print(", ρ=", ρ)
-            KA.@print(", q.tot=", q.tot)
-            KA.@print(", q.liq=", q.liq)
-            KA.@print(", q.ice=", q.ice)
-            KA.@print(", T=", sol.root)
-            KA.@print(", maxiter=", maxiter)
-            KA.@print(", tol=", tol.tol, "\n")
-        end
-        if error_on_non_convergence()
-            error("Failed to converge with printed set of inputs.")
-        end
+    if warn_msg(logger) && !sol.converged
+        KA.@print("-----------------------------------------\n")
+        KA.@print("maxiter reached in air_temperature_given_ρθq_nonlinear:\n")
+        KA.@print("    Method=SecantMethod")
+        KA.@print(", θ_liq_ice=", θ_liq_ice)
+        KA.@print(", ρ=", ρ)
+        KA.@print(", q.tot=", q.tot)
+        KA.@print(", q.liq=", q.liq)
+        KA.@print(", q.ice=", q.ice)
+        KA.@print(", T=", sol.root)
+        KA.@print(", maxiter=", maxiter)
+        KA.@print(", tol=", tol.tol, "\n")
+    end
+    if error_msg(logger) && !sol.converged
+        error("Failed to converge with printed set of inputs.")
     end
     return sol.root
 end
diff --git a/src/states.jl b/src/states.jl
index cbe54ce2..2b922684 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -298,7 +298,10 @@ and, optionally
     See the [`Thermodynamics`](@ref) for options.
  - `T_guess` initial guess for temperature in saturation adjustment
 """
+@inline PhaseEquil_ρeq(param_set::APS, args...) =
+    PhaseEquil_ρeq(WarnAndErrorLogger(), param_set::APS, args...)
 @inline function PhaseEquil_ρeq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     e_int::FT,
@@ -313,6 +316,7 @@ and, optionally
     phase_type = PhaseEquil{FT}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment(
+        logger,
         sat_adjust_method,
         param_set,
         e_int,
@@ -333,6 +337,7 @@ end
 # and relative_temperature_tol. maxiter and relative_temperature_tol
 # should be in sync with the PhaseEquil(...) constructor
 @inline function PhaseEquil_dev_only(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     e_int::FT,
@@ -342,6 +347,7 @@ end
     sat_adjust_method::Type{NM} = RS.NewtonsMethod,
 ) where {FT <: Real, NM}
     return PhaseEquil_ρeq(
+        logger,
         param_set,
         ρ,
         e_int,
@@ -353,7 +359,7 @@ end
 end
 
 """
-    PhaseEquil_ρθq(param_set, ρ, θ_liq_ice, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
+    PhaseEquil_ρθq([logger, ]param_set, ρ, θ_liq_ice, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
 
@@ -365,7 +371,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
  - `maxiter` maximum iterations for saturation adjustment
  - `T_guess` initial guess for temperature in saturation adjustment
 """
+@inline PhaseEquil_ρθq(param_set::APS, args...) =
+    PhaseEquil_ρθq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -380,6 +390,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
     phase_type = PhaseEquil{FT}
     tol = RS.RelativeSolutionTolerance(relative_temperature_tol)
     T = saturation_adjustment_given_ρθq(
+        logger,
         param_set,
         ρ,
         θ_liq_ice,
@@ -396,7 +407,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
 end
 
 """
-    PhaseEquil_ρTq(param_set, ρ, T, q_tot)
+    PhaseEquil_ρTq([logger, ]param_set, ρ, T, q_tot)
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -405,7 +416,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
  - `T` temperature
  - `q_tot` total specific humidity
 """
+@inline PhaseEquil_ρTq(param_set::APS, args...) =
+    PhaseEquil_ρTq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_ρTq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     T::FT,
@@ -419,7 +434,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 end
 
 """
-    PhaseEquil_pTq(param_set, p, T, q_tot)
+    PhaseEquil_pTq([logger, ]param_set, p, T, q_tot)
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -428,7 +443,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
  - `T` temperature
  - `q_tot` total specific humidity
 """
+@inline PhaseEquil_pTq(param_set::APS, args...) =
+    PhaseEquil_pTq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_pTq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     T::FT,
@@ -443,7 +462,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 end
 
 """
-    PhaseEquil_peq(param_set, p, e_int, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
+    PhaseEquil_peq([logger, ]param_set, p, e_int, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -453,7 +472,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
  - `q_tot` total specific humidity
  - `T_guess` initial guess for temperature in saturation adjustment
 """
+@inline PhaseEquil_peq(param_set::APS, args...) =
+    PhaseEquil_peq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_peq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     e_int::FT,
@@ -469,6 +492,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
     phase_type = PhaseEquil{FT}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment_given_peq(
+        logger,
         sat_adjust_method,
         param_set,
         p,
@@ -486,7 +510,7 @@ end
 
 
 """
-    PhaseEquil_phq(param_set, p, h, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
+    PhaseEquil_phq([logger, ]param_set, p, h, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -496,7 +520,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
  - `q_tot` total specific humidity
  - `T_guess` initial guess for temperature in saturation adjustment
 """
+@inline PhaseEquil_phq(param_set::APS, args...) =
+    PhaseEquil_phq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_phq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     h::FT,
@@ -512,6 +540,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
     phase_type = PhaseEquil{FT}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment_given_phq(
+        logger,
         sat_adjust_method,
         param_set,
         p,
@@ -529,7 +558,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 end
 
 """
-    PhaseEquil_ρpq(param_set, ρ, p, q_tot[, perform_sat_adjust=true, maxiter, sat_adjust_method, T_guess])
+    PhaseEquil_ρpq([logger, ]param_set, ρ, p, q_tot[, perform_sat_adjust=true, maxiter, sat_adjust_method, T_guess])
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -547,7 +576,11 @@ TODO: change input argument order: perform_sat_adjust is
       unique to this constructor, so it should be last.
       (breaking change)
 """
+@inline PhaseEquil_ρpq(param_set::APS, args...) =
+    PhaseEquil_ρpq(WarnAndErrorLogger(), param_set::APS, args...)
+
 @inline function PhaseEquil_ρpq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     p::FT,
@@ -562,6 +595,7 @@ TODO: change input argument order: perform_sat_adjust is
     phase_type = PhaseEquil{FT}
     if perform_sat_adjust
         T = saturation_adjustment_ρpq(
+            logger,
             sat_adjust_method,
             param_set,
             ρ,
@@ -584,7 +618,7 @@ end
 
 
 """
-    PhaseEquil_pθq(param_set, θ_liq_ice, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
+    PhaseEquil_pθq([logger, ]param_set, θ_liq_ice, q_tot[, maxiter, relative_temperature_tol, sat_adjust_method, T_guess])
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
 
@@ -597,7 +631,11 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
  - `sat_adjust_method` the numerical method to use.
  - `T_guess` initial guess for temperature in saturation adjustment
 """
+PhaseEquil_pθq(param_set::APS, args...) =
+    PhaseEquil_pθq(WarningLogger(), param_set, args...)
+
 @inline function PhaseEquil_pθq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     θ_liq_ice::FT,
@@ -613,6 +651,7 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from:
     phase_type = PhaseEquil{FT}
     q_tot_safe = clamp(q_tot, FT(0), FT(1))
     T = saturation_adjustment_given_pθq(
+        logger,
         sat_adjust_method,
         param_set,
         p,
@@ -702,7 +741,11 @@ and, optionally
  - `relative_temperature_tol` potential temperature for non-linear equation solve
  - `maxiter` maximum iterations for non-linear equation solve
 """
+PhaseNonEquil_ρθq(param_set::APS, args...) =
+    PhaseNonEquil_ρθq(WarningLogger(), param_set, args...)
+
 @inline function PhaseNonEquil_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -713,6 +756,7 @@ and, optionally
     phase_type = PhaseNonEquil{FT}
     tol = RS.RelativeSolutionTolerance(relative_temperature_tol)
     T = air_temperature_given_ρθq_nonlinear(
+        logger,
         param_set,
         ρ,
         θ_liq_ice,
diff --git a/test/relations.jl b/test/relations.jl
index 23a133db..286fafa1 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -357,6 +357,7 @@ end
     ρ = FT(1)
     phase_type = PhaseEquil
     @test TD.saturation_adjustment(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         internal_energy_sat(param_set, 300.0, ρ, q_tot, phase_type),
@@ -368,6 +369,7 @@ end
     ) ≈ 300.0
     @test abs(
         TD.saturation_adjustment(
+            TD.WarnAndErrorLogger(),
             RS.NewtonsMethod,
             param_set,
             internal_energy_sat(param_set, 300.0, ρ, q_tot, phase_type),
@@ -383,6 +385,7 @@ end
     ρ = FT(0.1)
     @test isapprox(
         TD.saturation_adjustment(
+            TD.WarnAndErrorLogger(),
             RS.SecantMethod,
             param_set,
             internal_energy_sat(param_set, 200.0, ρ, q_tot, phase_type),
@@ -397,6 +400,7 @@ end
     )
     @test abs(
         TD.saturation_adjustment(
+            TD.WarnAndErrorLogger(),
             RS.NewtonsMethod,
             param_set,
             internal_energy_sat(param_set, 200.0, ρ, q_tot, phase_type),
@@ -882,6 +886,7 @@ end
     (; q_tot, q_liq, q_ice, q_pt, RH, e_kin, e_pot) = profiles
 
     @test_throws ErrorException TD.saturation_adjustment.(
+        TD.WarnAndErrorLogger(),
         RS.NewtonsMethod,
         param_set,
         e_int,
@@ -893,6 +898,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         e_int,
@@ -904,6 +910,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_peq.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         p,
@@ -916,6 +923,7 @@ end
 
     T_virt = T # should not matter: testing for non-convergence
     @test_throws ErrorException TD.temperature_and_humidity_given_TᵥρRH.(
+        TD.WarnAndErrorLogger(),
         param_set,
         T_virt,
         ρ,
@@ -926,6 +934,7 @@ end
     )
 
     @test_throws ErrorException TD.air_temperature_given_ρθq_nonlinear.(
+        TD.WarnAndErrorLogger(),
         param_set,
         ρ,
         θ_liq_ice,
@@ -935,6 +944,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_ρθq.(
+        TD.WarnAndErrorLogger(),
         param_set,
         ρ,
         θ_liq_ice,
@@ -945,6 +955,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_pθq.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         p,
@@ -956,6 +967,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_pθq.(
+        TD.WarnAndErrorLogger(),
         RS.NewtonsMethodAD,
         param_set,
         p,
@@ -967,6 +979,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_ρpq.(
+        TD.WarnAndErrorLogger(),
         RS.NewtonsMethodAD,
         param_set,
         ρ,
@@ -1175,6 +1188,7 @@ end
         # Accurate but expensive `PhaseNonEquil_ρθq` constructor (Non-linear temperature from θ_liq_ice)
         T_non_linear =
             TD.air_temperature_given_ρθq_nonlinear.(
+                TD.WarnAndErrorLogger(),
                 param_set,
                 ρ,
                 θ_liq_ice,
@@ -1304,6 +1318,7 @@ end
 
         T_rec_qpt_rec =
             TD.temperature_and_humidity_given_TᵥρRH.(
+                TD.WarnAndErrorLogger(),
                 param_set,
                 T_virt,
                 ρ,
@@ -1707,16 +1722,14 @@ end
         )
 end
 
-TD.solution_type() = RS.VerboseSolution()
 @testset "Test data collection" begin
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
     param_set = TP.ThermodynamicsParameters(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     (; ρ, e_int, q_tot) = profiles
-    ts = PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot)
+    ts = PhaseEquil_ρeq.(TD.VerboseLogger(TD.WarnAndErrorLogger()), param_set, ρ, e_int, q_tot)
     data = TD.DataCollection.get_data()
     TD.DataCollection.print_summary(data)
     TD.DataCollection.reset_stats()
 end
-TD.solution_type() = RS.CompactSolution()