diff --git a/integration_tests/utils/Cases.jl b/integration_tests/utils/Cases.jl
index 58cfb323e..4bff19f8b 100644
--- a/integration_tests/utils/Cases.jl
+++ b/integration_tests/utils/Cases.jl
@@ -175,11 +175,9 @@ function initialize_profiles(self::CasesBase{SoaresCase}, Gr::Grid, GMV::GridMea
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = theta[k]
         GMV.T.values[k] = theta[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = theta[k]
     end
 
     set_bcs(GMV.H, Gr)
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.T, Gr)
     satadjust(GMV)
 end
@@ -272,9 +270,7 @@ function initialize_profiles(self::CasesBase{Nieuwstadt}, Gr::Grid, GMV::GridMea
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = theta[k]
         GMV.T.values[k] = theta[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = theta[k]
     end
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.H, Gr)
     set_bcs(GMV.T, Gr)
     satadjust(GMV)
@@ -392,9 +388,7 @@ function initialize_profiles(self::CasesBase{BomexCase}, Gr::Grid, GMV::GridMean
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = thetal[k]
         GMV.T.values[k] = thetal[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = thetal[k]
     end
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
@@ -527,10 +521,8 @@ function initialize_profiles(self::CasesBase{life_cycle_Tan2018}, Gr::Grid, GMV:
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = thetal[k]
         GMV.T.values[k] = thetal[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = thetal[k]
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
@@ -677,10 +669,8 @@ function initialize_profiles(self::CasesBase{Rico}, Gr::Grid, GMV::GridMeanVaria
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = thetal[k]
         GMV.T.values[k] = thetal[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = thetal[k]
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
@@ -852,12 +842,9 @@ function initialize_profiles(self::CasesBase{TRMM_LBA}, Gr::Grid, GMV::GridMeanV
         GMV.H.values[k] =
             thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
 
-        GMV.THL.values[k] =
-            thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
         theta_rho[k] = theta_rho_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], qv)
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
     satadjust(GMV)
@@ -1117,12 +1104,8 @@ function initialize_profiles(self::CasesBase{ARM_SGP}, Gr::Grid, GMV::GridMeanVa
         GMV.T.values[k] = Theta[k] * exner_c(Ref.p0_half[k])
         GMV.H.values[k] =
             thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
-
-        GMV.THL.values[k] =
-            thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
@@ -1266,11 +1249,7 @@ function initialize_profiles(self::CasesBase{GATE_III}, Gr::Grid, GMV::GridMeanV
 
         GMV.H.values[k] =
             thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
-
-        GMV.THL.values[k] =
-            thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], 0.0, 0.0, latent_heat(GMV.T.values[k]))
     end
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.T, Gr)
@@ -1452,7 +1431,6 @@ function initialize_profiles(self::CasesBase{DYCOMS_RF01}, Gr::Grid, GMV::GridMe
         # thermodynamic variable profile (either entropy or thetal)
         # (calculated based on T and ql profiles.
         # Here we use Rd, cp and L constants as defined in TurbulenceConvection)
-        GMV.THL.values[k] = t_to_thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], GMV.QL.values[k], qi)
         GMV.H.values[k] = t_to_thetali_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], GMV.QL.values[k], qi)
 
 
@@ -1472,7 +1450,6 @@ function initialize_profiles(self::CasesBase{DYCOMS_RF01}, Gr::Grid, GMV::GridMe
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.QL, Gr)
     set_bcs(GMV.H, Gr)
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.T, Gr)
     set_bcs(GMV.B, Gr)
     return
@@ -1602,10 +1579,8 @@ function initialize_profiles(self::CasesBase{GABLS}, Gr::Grid, GMV::GridMeanVari
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = thetal[k]
         GMV.T.values[k] = thetal[k] * exner_c(Ref.p0_half[k]) # No water content
-        GMV.THL.values[k] = thetal[k]
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.V, Gr)
     set_bcs(GMV.QT, Gr)
@@ -1703,10 +1678,8 @@ function initialize_profiles(self::CasesBase{SP}, Gr::Grid, GMV::GridMeanVariabl
     @inbounds for k in real_center_indicies(Gr)
         GMV.H.values[k] = thetal[k]
         GMV.T.values[k] = thetal[k] * exner_c(Ref.p0_half[k])
-        GMV.THL.values[k] = thetal[k]
     end
 
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.U, Gr)
     set_bcs(GMV.V, Gr)
     set_bcs(GMV.QT, Gr)
@@ -1849,12 +1822,10 @@ function initialize_profiles(self::CasesBase{DryBubble}, Gr::Grid, GMV::GridMean
     cinterior = real_center_indicies(Gr)
     z_half_in = off_arr(Gr.z_half[cinterior])
     thetali[cinterior] = pyinterp(z_half_in, z_in, thetali_in)
-    GMV.THL.values .= thetali
     GMV.H.values .= thetali
     @inbounds for k in real_center_indicies(Gr)
         GMV.QT.values[k] = 0.0
     end
-    set_bcs(GMV.THL, Gr)
     set_bcs(GMV.QT, Gr)
     set_bcs(GMV.H, Gr)
     satadjust(GMV)
diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl
index 274c68de6..39ae4911a 100644
--- a/src/EDMF_Environment.jl
+++ b/src/EDMF_Environment.jl
@@ -96,7 +96,6 @@ end
 
 function update_EnvVar(self::EnvironmentThermodynamics, k, EnvVar::EnvironmentVariables, T, H, qt, ql, rho)
     EnvVar.T.values[k] = T
-    EnvVar.THL.values[k] = H
     EnvVar.H.values[k] = H
     EnvVar.QT.values[k] = qt
     EnvVar.QL.values[k] = ql
@@ -151,7 +150,7 @@ function saturation_adjustment(self::EnvironmentThermodynamics, EnvVar::Environm
             k,
             EnvVar,
             EnvVar.T.values[k],
-            EnvVar.THL.values[k],
+            EnvVar.H.values[k],
             EnvVar.QT.values[k],
             EnvVar.QL.values[k],
             EnvVar.QT.values[k] - EnvVar.QL.values[k],
diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl
index e4ebd9aa1..f36ff2d4d 100644
--- a/src/EDMF_Updrafts.jl
+++ b/src/EDMF_Updrafts.jl
@@ -42,14 +42,12 @@ function initialize(self::UpdraftVariables, GMV::GridMeanVariables)
             self.QL.values[i, k] = GMV.QL.values[k]
             self.H.values[i, k] = GMV.H.values[k]
             self.T.values[i, k] = GMV.T.values[k]
-            self.THL.values[i, k] = GMV.THL.values[k]
         end
 
         self.Area.values[i, kc_surf] = self.updraft_fraction / self.n_updrafts
     end
 
     set_bcs(self.QT, self.Gr)
-    set_bcs(self.THL, self.Gr)
     set_bcs(self.H, self.Gr)
     set_means(self, GMV)
 
@@ -159,7 +157,7 @@ function initialize_DryBubble(self::UpdraftVariables, GMV::GridMeanVariables, Re
                 # self.T.values[i,k] = sa.T
             else
                 self.Area.values[i, k] = 0.0 #self.updraft_fraction/self.n_updrafts
-                self.H.values[i, k] = GMV.THL.values[k]
+                self.H.values[i, k] = GMV.H.values[k]
                 self.T.values[i, k] = GMV.T.values[k]
             end
         end
@@ -167,7 +165,6 @@ function initialize_DryBubble(self::UpdraftVariables, GMV::GridMeanVariables, Re
 
     set_bcs(self.QT, self.Gr)
     set_bcs(self.H, self.Gr)
-    set_bcs(self.THL, self.Gr)
     set_bcs(self.W, self.Gr)
     set_bcs(self.T, self.Gr)
 
@@ -222,7 +219,6 @@ function set_means(self::UpdraftVariables, GMV::GridMeanVariables)
                 self.T.bulkvalues[k] += self.Area.values[i, k] * self.T.values[i, k] / self.Area.bulkvalues[k]
                 self.RH.bulkvalues[k] += self.Area.values[i, k] * self.RH.values[i, k] / self.Area.bulkvalues[k]
                 self.B.bulkvalues[k] += self.Area.values[i, k] * self.B.values[i, k] / self.Area.bulkvalues[k]
-                self.THL.bulkvalues[k] += self.Area.values[i, k] * self.THL.values[i, k] / self.Area.bulkvalues[k]
                 self.W.bulkvalues[k] += (
                     (self.Area.values[i, k] + self.Area.values[i, k + 1]) * self.W.values[i, k] /
                     (self.Area.bulkvalues[k] + self.Area.bulkvalues[k + 1])
@@ -235,7 +231,6 @@ function set_means(self::UpdraftVariables, GMV::GridMeanVariables)
             self.H.bulkvalues[k] = GMV.H.values[k]
             self.RH.bulkvalues[k] = GMV.RH.values[k]
             self.T.bulkvalues[k] = GMV.T.values[k]
-            self.THL.bulkvalues[k] = GMV.THL.values[k]
             self.B.bulkvalues[k] = 0.0
             self.W.bulkvalues[k] = 0.0
         end
@@ -261,7 +256,6 @@ function set_new_with_values(self::UpdraftVariables)
             self.QT.new[i, k] = self.QT.values[i, k]
             self.QL.new[i, k] = self.QL.values[i, k]
             self.H.new[i, k] = self.H.values[i, k]
-            self.THL.new[i, k] = self.THL.values[i, k]
             self.T.new[i, k] = self.T.values[i, k]
             self.B.new[i, k] = self.B.values[i, k]
         end
@@ -282,7 +276,6 @@ function set_old_with_values(self::UpdraftVariables)
             self.QT.old[i, k] = self.QT.values[i, k]
             self.QL.old[i, k] = self.QL.values[i, k]
             self.H.old[i, k] = self.H.values[i, k]
-            self.THL.old[i, k] = self.THL.values[i, k]
             self.T.old[i, k] = self.T.values[i, k]
             self.B.old[i, k] = self.B.values[i, k]
         end
@@ -303,7 +296,6 @@ function set_values_with_new(self::UpdraftVariables)
             self.QT.values[i, k] = self.QT.new[i, k]
             self.QL.values[i, k] = self.QL.new[i, k]
             self.H.values[i, k] = self.H.new[i, k]
-            self.THL.values[i, k] = self.THL.new[i, k]
             self.T.values[i, k] = self.T.new[i, k]
             self.B.values[i, k] = self.B.new[i, k]
         end
@@ -324,7 +316,6 @@ function io(self::UpdraftVariables, Stats::NetCDFIO_Stats, Ref::ReferenceState)
         write_profile(Stats, "updraft_thetal", self.H.bulkvalues[cinterior])
     else
         write_profile(Stats, "updraft_s", self.H.bulkvalues[cinterior])
-        #write_profile(Stats, "updraft_thetal", self.THL.bulkvalues[cinterior])
     end
 
     write_profile(Stats, "updraft_temperature", self.T.bulkvalues[cinterior])
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index f4e30ef81..4d2acd882 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -413,7 +413,7 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
 
         QT_cut = cut(self.EnvVar.QT.values, grid, k)
         grad_qt = c∇(QT_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
-        THL_cut = cut(self.EnvVar.THL.values, grid, k)
+        THL_cut = cut(self.EnvVar.H.values, grid, k)
         grad_thl = c∇(THL_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0))
 
         # g/theta_ref
@@ -509,7 +509,7 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl
                     self.EnvVar.T.values[k],
                 ) * (
                     (1.0 + (eps_vi - 1.0) * self.EnvVar.QT.values[k]) * grad_thl +
-                    (eps_vi - 1.0) * self.EnvVar.THL.values[k] * grad_qt
+                    (eps_vi - 1.0) * self.EnvVar.H.values[k] * grad_qt
                 )
             )
         N = sqrt(fmax(g / thv * grad_th_eff, 0.0))
@@ -1160,7 +1160,6 @@ function zero_area_fraction_cleanup(self::EDMF_PrognosticTKE, GMV::GridMeanVaria
                 self.UpdVar.QT.values[i, k] = GMV.QT.values[k]
                 self.UpdVar.T.values[i, k] = GMV.T.values[k]
                 self.UpdVar.QL.values[i, k] = GMV.QL.values[k]
-                self.UpdVar.THL.values[i, k] = GMV.THL.values[k]
             end
         end
 
@@ -1171,7 +1170,6 @@ function zero_area_fraction_cleanup(self::EDMF_PrognosticTKE, GMV::GridMeanVaria
             self.EnvVar.QT.values[k] = GMV.QT.values[k]
             self.EnvVar.T.values[k] = GMV.T.values[k]
             self.EnvVar.QL.values[k] = GMV.QL.values[k]
-            self.EnvVar.THL.values[k] = GMV.THL.values[k]
         end
     end
 
@@ -1185,14 +1183,12 @@ function set_subdomain_bcs(self::EDMF_PrognosticTKE)
     set_bcs(self.UpdVar.W, grid)
     set_bcs(self.UpdVar.Area, grid)
     set_bcs(self.UpdVar.H, grid)
-    set_bcs(self.UpdVar.THL, grid)
     set_bcs(self.UpdVar.QT, grid)
     set_bcs(self.UpdVar.T, grid)
     set_bcs(self.UpdVar.B, grid)
 
     set_bcs(self.EnvVar.W, grid)
     set_bcs(self.EnvVar.H, grid)
-    set_bcs(self.EnvVar.THL, grid)
     set_bcs(self.EnvVar.T, grid)
     set_bcs(self.EnvVar.QL, grid)
     set_bcs(self.EnvVar.QT, grid)
@@ -1592,7 +1588,6 @@ function update_GMV_ED(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::C
         end
     end
     set_bcs(GMV.QT, grid)
-    set_bcs(GMV.THL, grid)
     set_bcs(GMV.H, grid)
     set_bcs(GMV.U, grid)
     set_bcs(GMV.V, grid)
@@ -1624,7 +1619,7 @@ function compute_tke_buoy(self::EDMF_PrognosticTKE, GMV::GridMeanVariables)
         cpm = cpm_c(qt_cloudy)
         grad_thl_minus = grad_thl_plus
         grad_qt_minus = grad_qt_plus
-        grad_thl_plus = (self.EnvVar.THL.values[k + 1] - self.EnvVar.THL.values[k]) * grid.dzi
+        grad_thl_plus = (self.EnvVar.H.values[k + 1] - self.EnvVar.H.values[k]) * grid.dzi
         grad_qt_plus = (self.EnvVar.QT.values[k + 1] - self.EnvVar.QT.values[k]) * grid.dzi
         prefactor = Rd * exner_c(ref_state.p0_half[k]) / ref_state.p0_half[k]
         d_alpha_thetal_dry = prefactor * (1.0 + (eps_vi - 1.0) * qt_dry)
@@ -1695,8 +1690,6 @@ function update_GMV_diagnostics(self::EDMF_PrognosticTKE, GMV::GridMeanVariables
 
         qv = GMV.QT.values[k] - GMV.QL.values[k]
 
-        GMV.THL.values[k] = t_to_thetali_c(p0_half[k], GMV.T.values[k], GMV.QT.values[k], GMV.QL.values[k], 0.0)
-
         GMV.B.values[k] = (
             self.UpdVar.Area.bulkvalues[k] * self.UpdVar.B.bulkvalues[k] +
             (1.0 - self.UpdVar.Area.bulkvalues[k]) * self.EnvVar.B.values[k]
@@ -2503,7 +2496,7 @@ function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, opti
     elseif option == "thetal_maxgrad"
 
         @inbounds for k in real_center_indicies(self.Gr)
-            ∇θ_liq = ∇_upwind(GMV.THL.values, self.Gr, k)
+            ∇θ_liq = ∇_upwind(GMV.H.values, self.Gr, k)
             if ∇θ_liq > ∇θ_liq_max
                 ∇θ_liq_max = ∇θ_liq
                 self.zi = self.Gr.z[k]
diff --git a/src/Variables.jl b/src/Variables.jl
index ed7c0756f..72e292f02 100644
--- a/src/Variables.jl
+++ b/src/Variables.jl
@@ -88,7 +88,6 @@ function update(self::GridMeanVariables, TS::TimeStepping)
     set_bcs(self.V, grid)
     set_bcs(self.H, grid)
     set_bcs(self.QT, grid)
-    set_bcs(self.THL, grid)
     set_bcs(self.TKE, grid)
 
     set_bcs(self.QTvar, grid)
@@ -184,7 +183,6 @@ function satadjust(self::GridMeanVariables)
         self.QL.values[k] = sa.ql
         self.T.values[k] = sa.T
         qv = qt - sa.ql
-        self.THL.values[k] = t_to_thetali_c(p0, sa.T, qt, sa.ql, 0.0)
         rho = rho_c(p0, sa.T, qt, qv)
         self.B.values[k] = buoyancy_c(self.Ref.rho0_half[k], rho)
         self.RH.values[k] = relative_humidity_c(self.Ref.p0_half[k], qt, qt - qv, 0.0, self.T.values[k])
diff --git a/src/types.jl b/src/types.jl
index c2627e2bd..3304bc4d6 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -268,7 +268,6 @@ mutable struct UpdraftVariables{A1}
     QL::UpdraftVariable
     RH::UpdraftVariable
     H::UpdraftVariable
-    THL::UpdraftVariable
     T::UpdraftVariable
     B::UpdraftVariable
     prognostic::Bool
@@ -289,7 +288,6 @@ mutable struct UpdraftVariables{A1}
         QL = UpdraftVariable(grid, nu, "half", "scalar", "ql", "kg/kg")
         RH = UpdraftVariable(grid, nu, "half", "scalar", "RH", "%")
         H = UpdraftVariable(grid, nu, "half", "scalar", "thetal", "K")
-        THL = UpdraftVariable(grid, nu, "half", "scalar", "thetal", "K")
         T = UpdraftVariable(grid, nu, "half", "scalar", "temperature", "K")
         B = UpdraftVariable(grid, nu, "half", "scalar", "buoyancy", "m^2/s^3")
 
@@ -315,7 +313,6 @@ mutable struct UpdraftVariables{A1}
             QL,
             RH,
             H,
-            THL,
             T,
             B,
             prognostic,
@@ -347,7 +344,6 @@ Base.@kwdef mutable struct GridMeanVariables{PS}
     QL::VariableDiagnostic
     T::VariableDiagnostic
     B::VariableDiagnostic
-    THL::VariableDiagnostic
     cloud_fraction::VariableDiagnostic
     EnvThermo_scheme::String
     TKE::VariableDiagnostic
@@ -379,7 +375,6 @@ function GridMeanVariables(namelist, Gr::Grid, Ref::ReferenceState, param_set::P
     QL = VariableDiagnostic(Gr, "half", "scalar", "sym", "ql", "kg/kg")
     T = VariableDiagnostic(Gr, "half", "scalar", "sym", "temperature", "K")
     B = VariableDiagnostic(Gr, "half", "scalar", "sym", "buoyancy", "m^2/s^3")
-    THL = VariableDiagnostic(Gr, "half", "scalar", "sym", "thetal", "K")
 
     cloud_fraction = VariableDiagnostic(Gr, "half", "scalar", "sym", "cloud fraction", "-")
 
@@ -413,7 +408,6 @@ function GridMeanVariables(namelist, Gr::Grid, Ref::ReferenceState, param_set::P
         QL,
         T,
         B,
-        THL,
         cloud_fraction,
         EnvThermo_scheme,
         TKE,
@@ -527,7 +521,6 @@ Base.@kwdef mutable struct EnvironmentVariables{PS}
     QT::EnvironmentVariable
     QL::EnvironmentVariable
     H::EnvironmentVariable
-    THL::EnvironmentVariable
     RH::EnvironmentVariable
     T::EnvironmentVariable
     B::EnvironmentVariable
@@ -548,7 +541,6 @@ function EnvironmentVariables(namelist, Gr::Grid, param_set::PS) where {PS}
     QL = EnvironmentVariable(Gr, "half", "scalar", "ql", "kg/kg")
     RH = EnvironmentVariable(Gr, "half", "scalar", "RH", "%")
     H = EnvironmentVariable(Gr, "half", "scalar", "thetal", "K")
-    THL = EnvironmentVariable(Gr, "half", "scalar", "thetal", "K")
     T = EnvironmentVariable(Gr, "half", "scalar", "temperature", "K")
     B = EnvironmentVariable(Gr, "half", "scalar", "buoyancy", "m^2/s^3")
     Area = EnvironmentVariable(Gr, "half", "scalar", "env_area", "-")
@@ -569,7 +561,6 @@ function EnvironmentVariables(namelist, Gr::Grid, param_set::PS) where {PS}
         QT,
         QL,
         H,
-        THL,
         RH,
         T,
         B,