Use Thermodynamics.jl for reference profiles

Project.toml

@@ -22,4 +22,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"

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"] = 3.9046069426558216e+00
 best_mse["updraft_area"] = 2.8331099674735910e+04

src/ReferenceState.jl

@@ -44,30 +44,31 @@ specific initialization function defined in Initialization.pyx
 """
 function initialize(self::ReferenceState, Gr::Grid, Stats::NetCDFIO_Stats)
 
-    self.sg = t_to_entropy_c(self.Pg, self.Tg, self.qtg, 0.0, 0.0)
+    FT = eltype(Gr)
+    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)
+
+    # 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_pTq(param_set, p_, T_g, q_tot_g)
+        R_m = TD.gas_constant_air(ts)
+        T = TD.air_temperature(ts)
+        return -FT(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(Gr)
     p_half = center_field(Gr)
 
     # Perform the integration
-    # TODO: replace with OrdinaryDiffEq
-
     z_span = (Gr.zmin, Gr.zmax)
     @show z_span
     prob = ODEProblem(rhs, logp, z_span)
@@ -107,21 +108,20 @@ function initialize(self::ReferenceState, Gr::Grid, Stats::NetCDFIO_Stats)
     qv_half = center_field(Gr)
 
     # Compute reference state thermodynamic profiles
-
     @inbounds for k in center_indicies(Gr)
-        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_pTq(param_set, p_half_[k], T_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(Gr)
-        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])
+        ts = TD.PhaseEquil_pTq(param_set, p_[k], T_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
 
     # Now do a sanity check to make sure that the Reference State entropy profile is uniform following

src/TurbulenceConvection.jl

@@ -9,6 +9,10 @@ import LambertW
 using Distributions: Normal, quantile
 using OffsetArrays
 using FastGaussQuadrature: gausshermite
+import Thermodynamics
+const TD = Thermodynamics
+using CLIMAParameters: AbstractEarthParameterSet
+using CLIMAParameters.Planet: grav
 using LinearAlgebra
 
 import CLIMAParameters