Merge #87

87: Remove THL, and replace with H r=charleskawczynski a=charleskawczynski Closes #59. Note that this may actually be behavior-changing if H and THL are ever out of sync. Co-authored-by: Charles Kawczynski <[email protected]>
CliMA · Aug 6, 2021 · 037be4f · 037be4f
2 parents 6ace5d2 + 656a1d6
commit 037be4f
Show file tree

Hide file tree

Showing 6 changed files with 6 additions and 63 deletions.
diff --git 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)
@@ -1118,12 +1105,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)
@@ -1267,11 +1250,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)
@@ -1453,7 +1432,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)
 
 
@@ -1473,7 +1451,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
@@ -1603,10 +1580,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)
@@ -1704,10 +1679,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)
@@ -1850,12 +1823,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
@@ -99,7 +99,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
@@ -154,7 +153,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
@@ -44,14 +44,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)
 
@@ -161,15 +159,14 @@ 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
     end
 
     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)
 
@@ -224,7 +221,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])
@@ -237,7 +233,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
@@ -263,7 +258,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
@@ -284,7 +278,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
@@ -305,7 +298,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
@@ -326,7 +318,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
@@ -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
@@ -92,7 +92,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)
@@ -188,7 +187,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])