diff --git a/Project.toml b/Project.toml
index d98f7a3aad..130a66c985 100644
--- a/Project.toml
+++ b/Project.toml
@@ -21,4 +21,5 @@ PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
diff --git a/integration_tests/TRMM_LBA.jl b/integration_tests/TRMM_LBA.jl
index e98149ede6..3653161aa8 100644
--- a/integration_tests/TRMM_LBA.jl
+++ b/integration_tests/TRMM_LBA.jl
@@ -10,6 +10,9 @@ include(joinpath("utils", "generate_namelist.jl"))
 include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
+# Note: temperatures in this case become extremely low.
+CLIMAParameters.Planet.T_freeze(::EarthParameterSet) = 100.0
+
 best_mse = OrderedDict()
 best_mse["qt_mean"] = 1.7430143247852277e+00
 best_mse["updraft_area"] = 3.2544587990881948e+04
diff --git a/src/ReferenceState.jl b/src/ReferenceState.jl
index dacc85ae8c..b1d75720b4 100644
--- a/src/ReferenceState.jl
+++ b/src/ReferenceState.jl
@@ -44,30 +44,33 @@ specific initialization function defined in Initialization.pyx
 """
 function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats)
 
-    self.sg = t_to_entropy_c(self.Pg, self.Tg, self.qtg, 0.0, 0.0)
+    FT = eltype(grid)
+    P_g = self.Pg
+    T_g = self.Tg
+    q_tot_g = self.qtg
+    param_set = parameter_set(self)
+    q_pt_g = TD.PhasePartition(q_tot_g)
+    ts_g = TD.PhaseEquil_pTq(param_set, P_g, T_g, q_tot_g)
+    θ_liq_ice_g = TD.liquid_ice_pottemp(ts_g)
+
+    # We are integrating the log pressure so need to take the log of the
+    # surface pressure
+    logp = log(P_g)
 
     # Form a right hand side for integrating the hydrostatic equation to
     # determine the reference pressure
-
-    function rhs(p, u, z)
-        ret = eos(exp(p), self.qtg, self.sg; t_to_prog = t_to_entropy_c, prog_to_t = eos_first_guess_entropy)
-        q_i = 0.0
-        q_l = ret.ql
-        T = ret.T
-        return -g / (Rd * T * (1.0 - self.qtg + eps_vi * (self.qtg - q_l - q_i)))
+    function rhs(logp, u, z)
+        p_ = exp(logp)
+        ts = TD.PhaseEquil_pθq(param_set, p_, θ_liq_ice_g, q_tot_g)
+        R_m = TD.gas_constant_air(ts)
+        T = TD.air_temperature(ts)
+        return -FT(CPP.grav(param_set)) / (T * R_m)
     end
 
-    # We are integrating the log pressure so need to take the log of the
-    # surface pressure
-    p0 = log(self.Pg)
-    logp = p0
-
     p = face_field(grid)
     p_half = center_field(grid)
 
     # Perform the integration
-    # TODO: replace with OrdinaryDiffEq
-
     z_span = (grid.zmin, grid.zmax)
     @show z_span
     prob = ODEProblem(rhs, logp, z_span)
@@ -109,35 +112,22 @@ function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats)
     qv_half = center_field(grid)
 
     # Compute reference state thermodynamic profiles
-
     @inbounds for k in center_indicies(grid)
-        ret = eos(p_half_[k], self.qtg, self.sg; t_to_prog = t_to_entropy_c, prog_to_t = eos_first_guess_entropy)
-        temperature_half[k] = ret.T
-        ql_half[k] = ret.ql
-        qv_half[k] = self.qtg - (ql_half[k] + qi_half[k])
-        alpha_half[k] = alpha_c(p_half_[k], temperature_half[k], self.qtg, qv_half[k])
+        ts = TD.PhaseEquil_pθq(param_set, p_half_[k], θ_liq_ice_g, q_tot_g)
+        temperature_half[k] = TD.air_temperature(ts)
+        ql_half[k] = TD.liquid_specific_humidity(ts)
+        qv_half[k] = TD.vapor_specific_humidity(ts)
+        alpha_half[k] = TD.specific_volume(ts)
     end
 
     @inbounds for k in face_indicies(grid)
-        ret = eos(p_[k], self.qtg, self.sg; t_to_prog = t_to_entropy_c, prog_to_t = eos_first_guess_entropy)
-        temperature[k] = ret.T
-        ql[k] = ret.ql
-        qv[k] = self.qtg - (ql[k] + qi[k])
-        alpha[k] = alpha_c(p_[k], temperature[k], self.qtg, qv[k])
-    end
-
-    # Now do a sanity check to make sure that the Reference State entropy profile is uniform following
-    # saturation adjustment
-    local s
-    @inbounds for k in center_indicies(grid)
-        s = t_to_entropy_c(p_half[k], temperature_half[k], self.qtg, ql_half[k], qi_half[k])
-        if abs(s - self.sg) / self.sg > 0.01
-            println("Error in reference profiles entropy not constant !")
-            println("Likely error in saturation adjustment")
-        end
+        ts = TD.PhaseEquil_pθq(param_set, p_[k], θ_liq_ice_g, q_tot_g)
+        temperature[k] = TD.air_temperature(ts)
+        ql[k] = TD.liquid_specific_humidity(ts)
+        qv[k] = TD.vapor_specific_humidity(ts)
+        alpha[k] = TD.specific_volume(ts)
     end
 
-
     self.alpha0_half = alpha_half
     self.alpha0 = alpha
     self.p0 = p_
diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl
index 834879f476..df3ef6e07f 100644
--- a/src/TurbulenceConvection.jl
+++ b/src/TurbulenceConvection.jl
@@ -7,11 +7,15 @@ using LinearAlgebra
 import Dierckx
 import Statistics
 import LambertW
+import Thermodynamics
 import Distributions
 import FastGaussQuadrature
-
 import CLIMAParameters
+
+const TD = Thermodynamics
+
 const CP = CLIMAParameters
+const CPP = CP.Planet
 const APS = CP.AbstractEarthParameterSet
 
 const CPMP = CP.Atmos.Microphysics