diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl
index f0be33d9ba..d6e9d71f80 100644
--- a/src/EDMF_Environment.jl
+++ b/src/EDMF_Environment.jl
@@ -145,7 +145,7 @@ function saturation_adjustment(self::EnvironmentThermodynamics, EnvVar::Environm
     mph = mph_struct()
 
     @inbounds for k in real_center_indicies(self.Gr)
-        sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
+        sa = eos(self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
 
         EnvVar.T.values[k] = sa.T
         EnvVar.QL.values[k] = sa.ql
@@ -184,7 +184,7 @@ function sgs_mean(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables,
 
     @inbounds for k in real_center_indicies(self.Gr)
         # condensation
-        sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
+        sa = eos(self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
         # autoconversion and accretion
         mph = microphysics_rain_src(
             Rain.rain_model,
@@ -325,7 +325,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari
                     end
 
                     # condensation
-                    sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, self.Ref.p0_half[k], qt_hat, h_hat)
+                    sa = eos(self.Ref.p0_half[k], qt_hat, h_hat)
                     # autoconversion and accretion
                     mph = microphysics_rain_src(
                         Rain.rain_model,
@@ -411,7 +411,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari
 
         else
             # if variance and covariance are zero do the same as in SA_mean
-            sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
+            sa = eos(self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k])
             mph = microphysics_rain_src(
                 Rain.rain_model,
                 Rain.max_supersaturation,
diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl
index 0b4736d001..ddca5715b8 100644
--- a/src/EDMF_Updrafts.jl
+++ b/src/EDMF_Updrafts.jl
@@ -153,8 +153,6 @@ function initialize_DryBubble(self::UpdraftVariables, GMV::GridMeanVariables, Re
                 # for now temperature is provided as diagnostics from LES
 
                 # sa = eos(
-                #     t_to_thetali_c,
-                #     eos_first_guess_thetal,
                 #     Ref.p0_half[k],
                 #     self.QT.values[i,k],
                 #     self.H.values[i,k]
@@ -439,7 +437,7 @@ function buoyancy(
                 elseif UpdVar.Area.values[i, k - 1] > 0.0 && k > self.Gr.gw
                     # TODO: report bug:
                     # qt and h were not defined here before the function call.
-                    sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, self.Ref.p0_half[k], qt, h)
+                    sa = eos(self.Ref.p0_half[k], qt, h)
                     qt -= sa.ql
                     qv = qt
                     t = sa.T
diff --git a/src/ReferenceState.jl b/src/ReferenceState.jl
index 1be4606fc2..97485e4823 100644
--- a/src/ReferenceState.jl
+++ b/src/ReferenceState.jl
@@ -50,7 +50,7 @@ function initialize(self::ReferenceState, Gr::Grid, Stats::NetCDFIO_Stats)
     # determine the reference pressure
 
     function rhs(p, u, z)
-        ret = eos(t_to_entropy_c, eos_first_guess_entropy, exp(p), self.qtg, self.sg)
+        ret = eos(exp(p), self.qtg, self.sg)
         q_i = 0.0
         q_l = ret.ql
         T = ret.T
@@ -109,7 +109,7 @@ function initialize(self::ReferenceState, Gr::Grid, Stats::NetCDFIO_Stats)
     # Compute reference state thermodynamic profiles
 
     @inbounds for k in center_indicies(Gr)
-        ret = eos(t_to_entropy_c, eos_first_guess_entropy, p_half_[k], self.qtg, self.sg)
+        ret = eos(p_half_[k], self.qtg, self.sg)
         temperature_half[k] = ret.T
         ql_half[k] = ret.ql
         qv_half[k] = self.qtg - (ql_half[k] + qi_half[k])
@@ -117,7 +117,7 @@ function initialize(self::ReferenceState, Gr::Grid, Stats::NetCDFIO_Stats)
     end
 
     @inbounds for k in face_indicies(Gr)
-        ret = eos(t_to_entropy_c, eos_first_guess_entropy, p_[k], self.qtg, self.sg)
+        ret = eos(p_[k], self.qtg, self.sg)
         temperature[k] = ret.T
         ql[k] = ret.ql
         qv[k] = self.qtg - (ql[k] + qi[k])
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index ab651903e1..32607c6582 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -1372,13 +1372,7 @@ function solve_updraft_scalars(self::EDMF_PrognosticTKE, GMV::GridMeanVariables)
                 end
 
                 # saturation adjustment
-                sa = eos(
-                    self.UpdThermo.t_to_prog_fp,
-                    self.UpdThermo.prog_to_t_fp,
-                    ref_state.p0_half[k],
-                    self.UpdVar.QT.new[i, k],
-                    self.UpdVar.H.new[i, k],
-                )
+                sa = eos(ref_state.p0_half[k], self.UpdVar.QT.new[i, k], self.UpdVar.H.new[i, k])
                 self.UpdVar.QL.new[i, k] = sa.ql
                 self.UpdVar.T.new[i, k] = sa.T
                 continue
@@ -1427,13 +1421,7 @@ function solve_updraft_scalars(self::EDMF_PrognosticTKE, GMV::GridMeanVariables)
             end
 
             # saturation adjustment
-            sa = eos(
-                self.UpdThermo.t_to_prog_fp,
-                self.UpdThermo.prog_to_t_fp,
-                ref_state.p0_half[k],
-                self.UpdVar.QT.new[i, k],
-                self.UpdVar.H.new[i, k],
-            )
+            sa = eos(ref_state.p0_half[k], self.UpdVar.QT.new[i, k], self.UpdVar.H.new[i, k])
             self.UpdVar.QL.new[i, k] = sa.ql
             self.UpdVar.T.new[i, k] = sa.T
         end
diff --git a/src/Variables.jl b/src/Variables.jl
index db92703889..90068537cf 100644
--- a/src/Variables.jl
+++ b/src/Variables.jl
@@ -192,7 +192,7 @@ function satadjust(self::GridMeanVariables)
         h = self.H.values[k]
         qt = self.QT.values[k]
         p0 = self.Ref.p0_half[k]
-        sa = eos(self.t_to_prog_fp, self.prog_to_t_fp, p0, qt, h)
+        sa = eos(p0, qt, h)
         self.QL.values[k] = sa.ql
         self.T.values[k] = sa.T
         qv = qt - sa.ql
diff --git a/src/thermodynamic_functions.jl b/src/thermodynamic_functions.jl
index a8d30dc02b..61d38e5539 100644
--- a/src/thermodynamic_functions.jl
+++ b/src/thermodynamic_functions.jl
@@ -107,7 +107,7 @@ function eos_first_guess_entropy(H, pd, pv, qt)
     )
 end
 
-function eos(t_to_prog, prog_to_t, p0, qt, prog)
+function eos(p0, qt, θ_li)
     qv = qt
     ql = 0.0
 
@@ -115,7 +115,7 @@ function eos(t_to_prog, prog_to_t, p0, qt, prog)
 
     pv_1 = pv_c(p0, qt, qt)
     pd_1 = p0 - pv_1
-    T_1 = prog_to_t(prog, pd_1, pv_1, qt)
+    T_1 = eos_first_guess_thetal(θ_li, pd_1, pv_1, qt)
     pv_star_1 = pv_star(T_1)
     qv_star_1 = qv_star_c(p0, qt, pv_star_1)
 
@@ -127,8 +127,8 @@ function eos(t_to_prog, prog_to_t, p0, qt, prog)
 
     else
         ql_1 = qt - qv_star_1
-        prog_1 = t_to_prog(p0, T_1, qt, ql_1, 0.0)
-        f_1 = prog - prog_1
+        θ_li_1 = t_to_thetali_c(p0, T_1, qt, ql_1, 0.0)
+        f_1 = θ_li - θ_li_1
         T_2 = T_1 + ql_1 * latent_heat(T_1) / ((1.0 - qt) * cpd + qv_star_1 * cpv)
         delta_T = fabs(T_2 - T_1)
         qv_star_2 = 0
@@ -138,8 +138,8 @@ function eos(t_to_prog, prog_to_t, p0, qt, prog)
             pv_2 = pv_c(p0, qt, qv_star_2)
             pd_2 = p0 - pv_2
             ql_2 = qt - qv_star_2
-            prog_2 = t_to_prog(p0, T_2, qt, ql_2, 0.0)
-            f_2 = prog - prog_2
+            θ_li_2 = t_to_thetali_c(p0, T_2, qt, ql_2, 0.0)
+            f_2 = θ_li - θ_li_2
             T_n = T_2 - f_2 * (T_2 - T_1) / (f_2 - f_1)
             T_1 = T_2
             T_2 = T_n
diff --git a/src/types.jl b/src/types.jl
index f65529e53c..0a14e5d031 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -386,8 +386,6 @@ Base.@kwdef mutable struct GridMeanVariables{PS}
     QT::VariablePrognostic
     RH::VariablePrognostic
     H::VariablePrognostic
-    t_to_prog_fp::Function
-    prog_to_t_fp::Function
     QL::VariableDiagnostic
     T::VariableDiagnostic
     B::VariableDiagnostic
@@ -420,8 +418,6 @@ function GridMeanVariables(namelist, Gr::Grid, Ref::ReferenceState, param_set::P
     RH = VariablePrognostic(Gr, "half", "scalar", "sym", "RH", "%")
 
     H = VariablePrognostic(Gr, "half", "scalar", "sym", "thetal", "K")
-    t_to_prog_fp = t_to_thetali_c
-    prog_to_t_fp = eos_first_guess_thetal
 
     # Diagnostic Variables--same class as the prognostic variables, but we append to diagnostics list
     QL = VariableDiagnostic(Gr, "half", "scalar", "sym", "ql", "kg/kg")
@@ -478,8 +474,6 @@ function GridMeanVariables(namelist, Gr::Grid, Ref::ReferenceState, param_set::P
         QT,
         RH,
         H,
-        t_to_prog_fp,
-        prog_to_t_fp,
         QL,
         T,
         B,
@@ -503,8 +497,6 @@ struct UpdraftThermodynamics{A1, A2}
     Gr::Grid
     Ref::ReferenceState
     n_updraft::Int
-    t_to_prog_fp::Function
-    prog_to_t_fp::Function
     prec_source_h::A2
     prec_source_qt::A2
     prec_source_h_tot::A1
@@ -516,8 +508,6 @@ struct UpdraftThermodynamics{A1, A2}
         UpdVar::UpdraftVariables,
         Rain::RainVariables,
     )
-        t_to_prog_fp = t_to_thetali_c
-        prog_to_t_fp = eos_first_guess_thetal
 
         # rain source from each updraft from all sub-timesteps
         prec_source_h = center_field(Gr, n_updraft)
@@ -528,17 +518,7 @@ struct UpdraftThermodynamics{A1, A2}
         prec_source_qt_tot = center_field(Gr)
         A1 = typeof(prec_source_h_tot)
         A2 = typeof(prec_source_h)
-        return new{A1, A2}(
-            Gr,
-            Ref,
-            n_updraft,
-            t_to_prog_fp,
-            prog_to_t_fp,
-            prec_source_h,
-            prec_source_qt,
-            prec_source_h_tot,
-            prec_source_qt_tot,
-        )
+        return new{A1, A2}(Gr, Ref, n_updraft, prec_source_h, prec_source_qt, prec_source_h_tot, prec_source_qt_tot)
     end
 
 end
@@ -708,8 +688,6 @@ struct EnvironmentThermodynamics{A1}
     Ref::ReferenceState
     quadrature_order::Int
     quadrature_type::String
-    t_to_prog_fp::Function
-    prog_to_t_fp::Function
     qt_dry::A1
     th_dry::A1
     t_cloudy::A1
@@ -738,8 +716,6 @@ struct EnvironmentThermodynamics{A1}
         catch
             "gaussian"
         end
-        t_to_prog_fp = t_to_thetali_c
-        prog_to_t_fp = eos_first_guess_thetal
 
         qt_dry = center_field(Gr)
         th_dry = center_field(Gr)
@@ -761,8 +737,6 @@ struct EnvironmentThermodynamics{A1}
             Ref,
             quadrature_order,
             quadrature_type,
-            t_to_prog_fp,
-            prog_to_t_fp,
             qt_dry,
             th_dry,
             t_cloudy,