diff --git a/Project.toml b/Project.toml
index 49b440cc..8583904b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.10.2"
+version = "0.11.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
diff --git a/src/Parameters.jl b/src/Parameters.jl
index 40d2f1f3..b063691d 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -21,6 +21,7 @@ param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...);
 Base.@kwdef struct ThermodynamicsParameters{FT}
     T_0::FT
     MSLP::FT
+    p_ref_theta::FT
     cp_v::FT
     cp_l::FT
     cp_i::FT
diff --git a/src/isentropic.jl b/src/isentropic.jl
index f49d57a0..7543f951 100644
--- a/src/isentropic.jl
+++ b/src/isentropic.jl
@@ -28,10 +28,10 @@ function air_pressure_given_θ(
     Φ::FT,
     ::DryAdiabaticProcess,
 ) where {FT <: AbstractFloat}
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     _R_d::FT = TP.R_d(param_set)
     _cp_d::FT = TP.cp_d(param_set)
-    return MSLP * (1 - Φ / (θ * _cp_d))^(_cp_d / _R_d)
+    return p0 * (1 - Φ / (θ * _cp_d))^(_cp_d / _R_d)
 end
 
 """
@@ -72,6 +72,6 @@ function air_temperature(
 ) where {FT <: AbstractFloat}
     _R_d::FT = TP.R_d(param_set)
     _cp_d::FT = TP.cp_d(param_set)
-    MSLP::FT = TP.MSLP(param_set)
-    return (p / MSLP)^(_R_d / _cp_d) * θ
+    p0::FT = TP.p_ref_theta(param_set)
+    return (p / p0)^(_R_d / _cp_d) * θ
 end
diff --git a/src/relations.jl b/src/relations.jl
index 38e0583d..9c6e76f3 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -2293,12 +2293,12 @@ function air_temperature_given_ρθq(
     q::PhasePartition{FT} = q_pt_0(FT),
 ) where {FT <: Real}
 
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     cvm = cv_m(param_set, q)
     cpm = cp_m(param_set, q)
     R_m = gas_constant_air(param_set, q)
     κ = 1 - cvm / cpm
-    T_u = (ρ * R_m * θ_liq_ice / MSLP)^(R_m / cvm) * θ_liq_ice
+    T_u = (ρ * R_m * θ_liq_ice / p0)^(R_m / cvm) * θ_liq_ice
     T_1 = latent_heat_liq_ice(param_set, q) / cvm
     T_2 = -κ / (2 * T_u) * (latent_heat_liq_ice(param_set, q) / cvm)^2
     return T_u + T_1 + T_2
@@ -2545,13 +2545,13 @@ function exner_given_pressure(
     p::FT,
     q::PhasePartition{FT} = q_pt_0(FT),
 ) where {FT <: Real}
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     # gas constant and isobaric specific heat of moist air
     _R_m = gas_constant_air(param_set, q)
     _cp_m = cp_m(param_set, q)
 
-    # return (p / MSLP)^(_R_m / _cp_m)
-    return pow_hack(p / MSLP, _R_m / _cp_m)
+    # return (p / p0)^(_R_m / _cp_m)
+    return pow_hack(p / p0, _R_m / _cp_m)
 end
 
 """
diff --git a/test/TemperatureProfiles.jl b/test/TemperatureProfiles.jl
index 8f212c24..96555fa6 100644
--- a/test/TemperatureProfiles.jl
+++ b/test/TemperatureProfiles.jl
@@ -28,7 +28,7 @@ parameter_set(::Type{Float32}) = param_set_Float32
         param_set = parameter_set(FT)
         _grav = FT(TP.grav(param_set))
         _R_d = FT(TP.R_d(param_set))
-        _MSLP = FT(TP.MSLP(param_set))
+        _p_ref = FT(TP.MSLP(param_set))
 
         z = collect(range(FT(0), stop = FT(25e3), length = 100))
         n = 7
@@ -61,19 +61,19 @@ parameter_set(::Type{Float32}) = param_set_Float32
                 p = last.(args)
 
                 # Test that surface pressure is equal to the
-                # specified boundary condition (MSLP)
+                # specified boundary condition p_ref (by default set to MSLP)
                 mask_z_0 = z .≈ 0
-                @test all(p[mask_z_0] .≈ _MSLP)
+                @test all(p[mask_z_0] .≈ _p_ref)
 
-                function log_p_over_MSLP(_z)
+                function log_p_over_p_ref(_z)
                     _, _p = profile(param_set, _z)
-                    return log(_p / _MSLP)
+                    return log(_p / _p_ref)
                 end
-                ∇log_p_over_MSLP =
-                    _z -> ForwardDiff.derivative(log_p_over_MSLP, _z)
+                ∇log_p_over_p_ref =
+                    _z -> ForwardDiff.derivative(log_p_over_p_ref, _z)
                 # Uses density computed from pressure derivative
                 # and ideal gas law to reconstruct virtual temperature
-                T_virt_rec = -_grav ./ (_R_d .* ∇log_p_over_MSLP.(z))
+                T_virt_rec = -_grav ./ (_R_d .* ∇log_p_over_p_ref.(z))
                 @test all(T_virt_rec .≈ T_virt)
             end
         end
diff --git a/test/relations.jl b/test/relations.jl
index 0680ab1b..fe615244 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -86,7 +86,7 @@ compare_moisture(param_set, ts::PhaseNonEquil, q_pt::PhasePartition) = all((
         _T_triple = FT(TP.T_triple(param_set))
         _T_freeze = FT(TP.T_freeze(param_set))
         _T_min = FT(TP.T_min(param_set))
-        _MSLP = FT(TP.MSLP(param_set))
+        _p_ref_theta = FT(TP.p_ref_theta(param_set))
         _T_max = FT(TP.T_max(param_set))
         _kappa_d = FT(TP.kappa_d(param_set))
 
@@ -105,13 +105,13 @@ compare_moisture(param_set, ts::PhaseNonEquil, q_pt::PhasePartition) = all((
         # TODO: Use reasonable values for ambient temperature/pressure
         T∞, p∞ = T .* perturbation, p .* perturbation
         @test air_temperature.(param_set, p, θ_liq_ice, DryAdiabaticProcess()) ≈
-              (p ./ _MSLP) .^ (_R_d / _cp_d) .* θ_liq_ice
+              (p ./ _p_ref_theta) .^ (_R_d / _cp_d) .* θ_liq_ice
         @test TD.air_pressure_given_θ.(
             param_set,
             θ_liq_ice,
             Φ,
             DryAdiabaticProcess(),
-        ) ≈ _MSLP .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d)
+        ) ≈ _p_ref_theta .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d)
         @test air_pressure.(param_set, T, T∞, p∞, DryAdiabaticProcess()) ≈
               p∞ .* (T ./ T∞) .^ (FT(1) / _kappa_d)
     end
@@ -142,7 +142,7 @@ end
     _T_triple = FT(TP.T_triple(param_set))
     _T_freeze = FT(TP.T_freeze(param_set))
     _T_min = FT(TP.T_min(param_set))
-    _MSLP = FT(TP.MSLP(param_set))
+    _p_ref_theta = FT(TP.p_ref_theta(param_set))
     _T_max = FT(TP.T_max(param_set))
     _kappa_d = FT(TP.kappa_d(param_set))
     _T_icenuc = FT(TP.T_icenuc(param_set))
@@ -446,13 +446,16 @@ end
 
     # potential temperatures
     T = FT(300)
-    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _MSLP) === T
-    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _MSLP / 10) ≈
-          T * 10^(_R_d / _cp_d)
+    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _p_ref_theta) === T
     @test TD.liquid_ice_pottemp_given_pressure(
         param_set,
         T,
-        _MSLP / 10,
+        _p_ref_theta / 10,
+    ) ≈ T * 10^(_R_d / _cp_d)
+    @test TD.liquid_ice_pottemp_given_pressure(
+        param_set,
+        T,
+        _p_ref_theta / 10,
         PhasePartition(FT(1)),
     ) ≈ T * 10^(_R_v / _cp_v)
 
@@ -833,7 +836,7 @@ end
         @test count(air_temperature.(param_set, ts_upper) .== Ref(_T_freeze)) ≥
               217
         @test count(air_temperature.(param_set, ts_mid) .== Ref(_T_freeze)) ≥
-              1373
+              1296
         # we should do this instead, but we're failing because some inputs are bad
         # E.g. p ~ 110_000 Pa, q_tot ~ 0.16, which results in negative θ_liq_ice
         # This means that we should probably update our tested profiles.
@@ -992,7 +995,6 @@ end
     for ArrayType in array_types
         FT = eltype(ArrayType)
         param_set = parameter_set(FT)
-        _MSLP = FT(TP.MSLP(param_set))
 
         profiles = TestedProfiles.PhaseDryProfiles(param_set, ArrayType)
         @unpack T, p, e_int, h, ρ, θ_liq_ice, phase_type = profiles