diff --git a/docs/src/saturation_adjustment.jl b/docs/src/saturation_adjustment.jl
index d62481a6..4442e492 100644
--- a/docs/src/saturation_adjustment.jl
+++ b/docs/src/saturation_adjustment.jl
@@ -7,7 +7,6 @@ import Thermodynamics as TD
 const TP = TD.Parameters
 import CLIMAParameters as CP
 
-TD.print_warning() = false
 toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
 aliases = string.(fieldnames(TP.ThermodynamicsParameters))
 param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
@@ -59,8 +58,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,
@@ -68,10 +67,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 b47d7cf0..507521f8 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 = get_parameter_set(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..e4cf0a04 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(), )..., 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 63855122..5b75d740 100644
--- a/src/Thermodynamics.jl
+++ b/src/Thermodynamics.jl
@@ -58,26 +58,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.
-error_on_non_convergence() = true
-
-# Allow users to skip printing warnings on non-convergence
-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 7d4c214c..79cfd1a8 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -1491,6 +1491,7 @@ end
 
 """
     saturation_adjustment(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         e_int,
@@ -1524,6 +1525,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     e_int::FT,
@@ -1601,34 +1603,33 @@ function saturation_adjustment(
             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,
@@ -1664,6 +1665,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment_given_peq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -1706,26 +1708,24 @@ function saturation_adjustment_given_peq(
             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
@@ -1733,6 +1733,7 @@ end
 
 """
     saturation_adjustment_given_phq(
+        logger::ATL,
         sat_adjust_method,
         param_set,
         h,
@@ -1768,6 +1769,7 @@ using the given numerical method `sat_adjust_method`.
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment_given_phq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -1818,32 +1820,31 @@ function saturation_adjustment_given_phq(
             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,
         ρ,
@@ -1879,6 +1880,7 @@ using Newtons method using ForwardDiff.
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment_ρpq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     ρ::FT,
@@ -1914,26 +1916,24 @@ function saturation_adjustment_ρpq(
             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
@@ -1965,6 +1965,7 @@ end
 
 """
     saturation_adjustment_given_ρθq(
+        logger::ATL,
         param_set,
         ρ,
         θ_liq_ice,
@@ -1996,6 +1997,7 @@ by finding the root of
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment_given_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -2025,31 +2027,30 @@ function saturation_adjustment_given_ρθq(
                 q_tot,
             ) - θ_liq_ice,
         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,
@@ -2080,6 +2081,7 @@ by finding the root of
 See also [`saturation_adjustment`](@ref).
 """
 function saturation_adjustment_given_pθq(
+    logger::ATL,
     ::Type{sat_adjust_method},
     param_set::APS,
     p::FT,
@@ -2140,26 +2142,24 @@ function saturation_adjustment_given_pθq(
             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
@@ -2311,7 +2311,7 @@ function virt_temp_from_RH(
     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
 
@@ -2322,6 +2322,7 @@ The air temperature and `q_tot` where
  - `phase_type` a thermodynamic state type
 """
 function temperature_and_humidity_given_TᵥρRH(
+    logger::ATL,
     param_set::APS,
     T_virt::FT,
     ρ::FT,
@@ -2339,27 +2340,23 @@ function temperature_and_humidity_given_TᵥρRH(
             T_virt -
             virt_temp_from_RH(param_set, heavisided(T), ρ, RH, phase_type),
         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
 
@@ -2400,7 +2397,7 @@ function air_temperature_given_ρθq(
 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
 
@@ -2421,6 +2418,7 @@ by finding the root of
                                q) = 0`
 """
 function air_temperature_given_ρθq_nonlinear(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -2439,29 +2437,25 @@ function air_temperature_given_ρθq_nonlinear(
                 q,
             ),
         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 cd901186..9bcefeea 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -280,7 +280,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...)
 function PhaseEquil_ρeq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     e_int::FT,
@@ -296,6 +299,7 @@ function PhaseEquil_ρeq(
     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,
@@ -316,6 +320,7 @@ end
 # and relative_temperature_tol. maxiter and relative_temperature_tol
 # should be in sync with the PhaseEquil(...) constructor
 function PhaseEquil_dev_only(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     e_int::FT,
@@ -325,6 +330,7 @@ function PhaseEquil_dev_only(
     sat_adjust_method::Type{NM} = RS.NewtonsMethod,
 ) where {FT <: Real, NM}
     return PhaseEquil_ρeq(
+        logger,
         param_set,
         ρ,
         e_int,
@@ -336,7 +342,7 @@ function PhaseEquil_dev_only(
 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:
 
@@ -348,7 +354,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...)
+
 function PhaseEquil_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -363,6 +373,7 @@ function PhaseEquil_ρθq(
     phase_type = PhaseEquil{FT}
     tol = RS.RelativeSolutionTolerance(relative_temperature_tol)
     T = saturation_adjustment_given_ρθq(
+        logger,
         param_set,
         ρ,
         θ_liq_ice,
@@ -379,7 +390,7 @@ function PhaseEquil_ρθq(
 end
 
 """
-    PhaseEquil_ρTq(param_set, ρ, T, q_tot)
+    PhaseEquil_ρTq([logger, ]param_set, ρ, T, q_tot)
 
 Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.
 
@@ -388,7 +399,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...)
+
 function PhaseEquil_ρTq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     T::FT,
@@ -402,7 +417,7 @@ function PhaseEquil_ρTq(
 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.
 
@@ -411,7 +426,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...)
+
 function PhaseEquil_pTq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     T::FT,
@@ -426,7 +445,7 @@ function PhaseEquil_pTq(
 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.
 
@@ -436,7 +455,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...)
+
 function PhaseEquil_peq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     e_int::FT,
@@ -452,6 +475,7 @@ function PhaseEquil_peq(
     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,
@@ -469,7 +493,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.
 
@@ -479,7 +503,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...)
+
 function PhaseEquil_phq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     h::FT,
@@ -495,6 +523,7 @@ function PhaseEquil_phq(
     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,
@@ -512,7 +541,7 @@ function PhaseEquil_phq(
 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.
 
@@ -530,7 +559,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...)
+
 function PhaseEquil_ρpq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     p::FT,
@@ -545,6 +578,7 @@ function PhaseEquil_ρpq(
     phase_type = PhaseEquil{FT}
     if perform_sat_adjust
         T = saturation_adjustment_ρpq(
+            logger,
             sat_adjust_method,
             param_set,
             ρ,
@@ -567,7 +601,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:
 
@@ -580,7 +614,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...)
+
 function PhaseEquil_pθq(
+    logger::ATL,
     param_set::APS,
     p::FT,
     θ_liq_ice::FT,
@@ -596,6 +634,7 @@ function PhaseEquil_pθq(
     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,
@@ -685,7 +724,10 @@ 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...)
 function PhaseNonEquil_ρθq(
+    logger::ATL,
     param_set::APS,
     ρ::FT,
     θ_liq_ice::FT,
@@ -696,6 +738,7 @@ function PhaseNonEquil_ρθq(
     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 fe615244..77deb5f2 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -366,6 +366,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),
@@ -377,6 +378,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),
@@ -392,6 +394,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),
@@ -406,6 +409,7 @@ end
     )
     @test abs(
         TD.saturation_adjustment(
+            TD.WarnAndErrorLogger(),
             RS.NewtonsMethod,
             param_set,
             internal_energy_sat(param_set, 200.0, ρ, q_tot, phase_type),
@@ -891,6 +895,7 @@ end
     @unpack 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,
@@ -902,6 +907,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         e_int,
@@ -913,6 +919,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_peq.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         p,
@@ -925,6 +932,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,
         ρ,
@@ -935,6 +943,7 @@ end
     )
 
     @test_throws ErrorException TD.air_temperature_given_ρθq_nonlinear.(
+        TD.WarnAndErrorLogger(),
         param_set,
         ρ,
         θ_liq_ice,
@@ -944,6 +953,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_ρθq.(
+        TD.WarnAndErrorLogger(),
         param_set,
         ρ,
         θ_liq_ice,
@@ -954,6 +964,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_pθq.(
+        TD.WarnAndErrorLogger(),
         RS.SecantMethod,
         param_set,
         p,
@@ -965,6 +976,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_given_pθq.(
+        TD.WarnAndErrorLogger(),
         RS.NewtonsMethodAD,
         param_set,
         p,
@@ -976,6 +988,7 @@ end
     )
 
     @test_throws ErrorException TD.saturation_adjustment_ρpq.(
+        TD.WarnAndErrorLogger(),
         RS.NewtonsMethodAD,
         param_set,
         ρ,
@@ -1184,6 +1197,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,
@@ -1313,6 +1327,7 @@ end
 
         T_rec_qpt_rec =
             TD.temperature_and_humidity_given_TᵥρRH.(
+                TD.WarnAndErrorLogger(),
                 param_set,
                 T_virt,
                 ρ,
@@ -1716,16 +1731,21 @@ end
         )
 end
 
-TD.solution_type() = RS.VerboseSolution()
 @testset "Test data collection" begin
     ArrayType = Array{Float64}
     FT = eltype(ArrayType)
     param_set = parameter_set(FT)
     profiles = TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
     @unpack ρ, 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()