Merge #25

25: Use Thermodynamics.jl for reference state profiles r=charleskawczynski a=charleskawczynski

This PR adds the use of Thermodynamics.jl to help construct the reference state profiles.

Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>

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"

docs/src/ReferenceStates.md

@@ -62,16 +62,19 @@ for case_name in (
     "life_cycle_Tan2018",
     "Soares",
     "Rico",
-    "TRMM_LBA",
     "ARM_SGP",
-    "GATE_III",
     "DYCOMS_RF01",
     "GABLS",
     "SP",
     "DryBubble",
 )
     export_ref_profile(case_name)
 end;
+
+# Note: temperatures in this case become extremely low.
+CLIMAParameters.Planet.T_freeze(::EarthParameterSet) = 100.0
+export_ref_profile("TRMM_LBA")
+export_ref_profile("GATE_III")
 ```
 
 ## Bomex

integration_tests/ARM_SGP.jl

@@ -11,18 +11,18 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 3.4761388491700224e-01
-best_mse["updraft_area"] = 1.9655632954763928e+03
-best_mse["updraft_w"] = 3.2377093752385053e+02
-best_mse["updraft_qt"] = 1.2452570404535185e+01
-best_mse["updraft_thetal"] = 2.7685659900036992e+01
-best_mse["u_mean"] = 8.7998547277817906e+01
-best_mse["tke_mean"] = 5.9501141635391093e+02
-best_mse["temperature_mean"] = 1.4210607813065063e-04
-best_mse["ql_mean"] = 2.1890556031513017e+02
-best_mse["thetal_mean"] = 1.3953340654098429e-04
-best_mse["Hvar_mean"] = 1.1148198912643536e+03
-best_mse["QTvar_mean"] = 3.0703058500395360e+02
+best_mse["qt_mean"] = 3.4297134939930601e-01
+best_mse["updraft_area"] = 1.9829008493073268e+03
+best_mse["updraft_w"] = 3.1072236638537851e+02
+best_mse["updraft_qt"] = 1.2698381273033208e+01
+best_mse["updraft_thetal"] = 2.7686390354792074e+01
+best_mse["u_mean"] = 8.7998547277817863e+01
+best_mse["tke_mean"] = 6.7319699431828644e+02
+best_mse["temperature_mean"] = 1.3960321509272764e-04
+best_mse["ql_mean"] = 2.9854035802671558e+02
+best_mse["thetal_mean"] = 1.4200441599983529e-04
+best_mse["Hvar_mean"] = 1.3755159256021079e+03
+best_mse["QTvar_mean"] = 2.6566877259844244e+02
 
 @testset "ARM_SGP" begin
     println("Running ARM_SGP...")

integration_tests/Bomex.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 1.1068396401839847e-01
-best_mse["updraft_area"] = 7.2515940450436165e+02
-best_mse["updraft_w"] = 2.5856202193762641e+01
-best_mse["updraft_qt"] = 4.0061911576818243e+00
-best_mse["updraft_thetal"] = 2.1549386071253458e+01
-best_mse["v_mean"] = 6.5501952940732622e+01
-best_mse["u_mean"] = 5.3295928431425537e+01
-best_mse["tke_mean"] = 3.9582886707067345e+01
-best_mse["temperature_mean"] = 4.5842727523070819e-05
-best_mse["ql_mean"] = 6.5924727808279036e+00
-best_mse["thetal_mean"] = 4.6222694206747763e-05
-best_mse["Hvar_mean"] = 2.5900864672158934e+02
-best_mse["QTvar_mean"] = 4.6280994372607324e+01
+best_mse["qt_mean"] = 9.5080504754167433e-02
+best_mse["updraft_area"] = 7.2534663873051477e+02
+best_mse["updraft_w"] = 2.6358909128208378e+01
+best_mse["updraft_qt"] = 4.0334302720334296e+00
+best_mse["updraft_thetal"] = 2.1547995304586468e+01
+best_mse["v_mean"] = 6.5231550206742128e+01
+best_mse["u_mean"] = 5.3292889729286458e+01
+best_mse["tke_mean"] = 3.8860563196392924e+01
+best_mse["temperature_mean"] = 3.8436817065445734e-05
+best_mse["ql_mean"] = 3.7878522086312807e+00
+best_mse["thetal_mean"] = 3.9063754163395865e-05
+best_mse["Hvar_mean"] = 7.2639655951863034e+01
+best_mse["QTvar_mean"] = 2.2435640414468629e+01
 
 @testset "Bomex" begin
     println("Running Bomex...")

integration_tests/DYCOMS_RF01.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 1.5987989288888820e-02
-best_mse["ql_mean"] = 5.9262500323990217e+00
-best_mse["updraft_area"] = 2.3887443639936907e+02
-best_mse["updraft_w"] = 4.3317426538007693e+00
-best_mse["updraft_qt"] = 1.1983341214387264e+00
-best_mse["updraft_thetal"] = 1.2739627177179765e+01
-best_mse["v_mean"] = 3.9736878190256498e+01
-best_mse["u_mean"] = 3.7037227839730683e+01
-best_mse["tke_mean"] = 1.4955295527017412e+01
-best_mse["temperature_mean"] = 1.6122517253407575e-05
-best_mse["thetal_mean"] = 1.7630642054888154e-05
-best_mse["Hvar_mean"] = 1.2452017842800817e+04
-best_mse["QTvar_mean"] = 7.6626218978684426e+02
+best_mse["qt_mean"] = 1.6029892735513345e-02
+best_mse["ql_mean"] = 6.0295296107176544e+00
+best_mse["updraft_area"] = 2.3847762474738320e+02
+best_mse["updraft_w"] = 4.3312794675046318e+00
+best_mse["updraft_qt"] = 1.1929325513431135e+00
+best_mse["updraft_thetal"] = 1.2740562957739254e+01
+best_mse["v_mean"] = 3.9737745703044638e+01
+best_mse["u_mean"] = 3.7038099069367689e+01
+best_mse["tke_mean"] = 1.4950275863942871e+01
+best_mse["temperature_mean"] = 1.6731316538438984e-05
+best_mse["thetal_mean"] = 1.7840937307304572e-05
+best_mse["Hvar_mean"] = 1.2023987118420042e+04
+best_mse["QTvar_mean"] = 7.5103635873377698e+02
 
 @testset "DYCOMS_RF01" begin
     println("Running DYCOMS_RF01...")

integration_tests/GABLS.jl

@@ -12,14 +12,14 @@ using .NameList
 
 best_mse = OrderedDict()
 
-best_mse["updraft_thetal"] = 5.0248682056386063e+00
-best_mse["v_mean"] = 4.4635954562618867e+00
-best_mse["u_mean"] = 9.6412415987917157e+00
-best_mse["tke_mean"] = 2.4692059869443228e+00
-best_mse["temperature_mean"] = 8.8268614667882937e-06
-best_mse["thetal_mean"] = 8.7932314577535832e-06
-best_mse["Hvar_mean"] = 1.2890693335381430e+01
-best_mse["QTvar_mean"] = 4.4413913923170839e-01
+best_mse["updraft_thetal"] = 5.0248694008424222e+00
+best_mse["v_mean"] = 4.4630163253493214e+00
+best_mse["u_mean"] = 9.6412500983148846e+00
+best_mse["tke_mean"] = 2.4683562232142533e+00
+best_mse["temperature_mean"] = 8.8678308948818863e-06
+best_mse["thetal_mean"] = 8.7900622849842488e-06
+best_mse["Hvar_mean"] = 1.2891667586703637e+01
+best_mse["QTvar_mean"] = 4.4424755799649029e-01
 
 @testset "GABLS" begin
     println("Running GABLS...")

integration_tests/GATE_III.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
+
 @testset "GATE_III" begin
     println("Running GATE_III...")
     namelist = default_namelist("GATE_III")

integration_tests/Nieuwstadt.jl

@@ -11,14 +11,14 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["updraft_area"] = 5.9103576007452034e+02
-best_mse["updraft_w"] = 2.6741082264407293e+01
-best_mse["updraft_thetal"] = 3.0475272117346940e+01
-best_mse["u_mean"] = 1.5245873874565478e+02
-best_mse["tke_mean"] = 7.3707728452173882e+01
-best_mse["temperature_mean"] = 1.2199914281992930e-05
-best_mse["thetal_mean"] = 1.2067109520894632e-05
-best_mse["Hvar_mean"] = 1.9140690462510639e+02
+best_mse["updraft_area"] = 5.9412545314262150e+02
+best_mse["updraft_w"] = 2.6952552348430249e+01
+best_mse["updraft_thetal"] = 3.0475271481866773e+01
+best_mse["u_mean"] = 1.5247657003047263e+02
+best_mse["tke_mean"] = 7.3600391357561151e+01
+best_mse["temperature_mean"] = 1.1971868506884125e-05
+best_mse["thetal_mean"] = 1.2117924470979204e-05
+best_mse["Hvar_mean"] = 1.8622261767371862e+02
 
 @testset "Nieuwstadt" begin
     println("Running Nieuwstadt...")

integration_tests/Rico.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 4.2359420114784679e-01
-best_mse["updraft_area"] = 1.6722742675049749e+03
-best_mse["updraft_w"] = 2.8300484466673169e+02
-best_mse["updraft_qt"] = 1.7904453430368505e+01
-best_mse["updraft_thetal"] = 6.3573804437983789e+01
-best_mse["v_mean"] = 1.0634096403180666e+02
-best_mse["u_mean"] = 1.1434656097951681e+02
-best_mse["tke_mean"] = 8.6824438228445820e+02
-best_mse["temperature_mean"] = 1.2576042829475203e-04
-best_mse["ql_mean"] = 8.9402683009257942e+01
-best_mse["thetal_mean"] = 1.1896802334339157e-04
-best_mse["Hvar_mean"] = 2.0228054123757493e+03
-best_mse["QTvar_mean"] = 8.7748179858993274e+02
+best_mse["qt_mean"] = 4.2888329142054360e-01
+best_mse["updraft_area"] = 1.5804026921690113e+03
+best_mse["updraft_w"] = 3.0412249337738064e+02
+best_mse["updraft_qt"] = 1.8041619418228048e+01
+best_mse["updraft_thetal"] = 6.3565634303932690e+01
+best_mse["v_mean"] = 1.0637419570945200e+02
+best_mse["u_mean"] = 1.1444385357013196e+02
+best_mse["tke_mean"] = 8.6459012010773438e+02
+best_mse["temperature_mean"] = 1.2577596756983134e-04
+best_mse["ql_mean"] = 5.0938207504953844e+01
+best_mse["thetal_mean"] = 1.2604124901901089e-04
+best_mse["Hvar_mean"] = 1.4598483199237567e+03
+best_mse["QTvar_mean"] = 7.0850350571186209e+02
 
 @testset "Rico" begin
     println("Running Rico...")

integration_tests/Soares.jl

@@ -11,16 +11,16 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 1.9338238625419904e-01
-best_mse["updraft_area"] = 4.8852467389339699e+02
-best_mse["updraft_w"] = 3.1250501928635639e+01
-best_mse["updraft_qt"] = 1.1136647539128006e+01
-best_mse["updraft_thetal"] = 2.2395194054335001e+01
-best_mse["u_mean"] = 7.2109940121868794e+02
-best_mse["tke_mean"] = 4.7919324815294324e+01
-best_mse["temperature_mean"] = 1.4374232147002715e-05
-best_mse["thetal_mean"] = 1.4379965874656967e-05
-best_mse["Hvar_mean"] = 2.3305425470797826e+02
+best_mse["qt_mean"] = 1.9368730351926020e-01
+best_mse["updraft_area"] = 4.8515130524256091e+02
+best_mse["updraft_w"] = 3.1260485555848359e+01
+best_mse["updraft_qt"] = 1.1137263871861222e+01
+best_mse["updraft_thetal"] = 2.2395183773350311e+01
+best_mse["u_mean"] = 7.2111916707127898e+02
+best_mse["tke_mean"] = 4.7882808390238374e+01
+best_mse["temperature_mean"] = 1.3818642271008120e-05
+best_mse["thetal_mean"] = 1.4383770217873919e-05
+best_mse["Hvar_mean"] = 2.3235278609470686e+02
 
 @testset "Soares" begin
     println("Running Soares...")

integration_tests/TRMM_LBA.jl

@@ -10,20 +10,23 @@ 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.7430285968163504e+00
-best_mse["updraft_area"] = 3.2545401197014067e+04
-best_mse["updraft_w"] = 9.2113782722003907e+02
-best_mse["updraft_qt"] = 3.0306781659879494e+01
-best_mse["updraft_thetal"] = 1.1108867837608678e+02
-best_mse["v_mean"] = 2.9291577336310036e+02
-best_mse["u_mean"] = 1.6874237352911189e+03
-best_mse["tke_mean"] = 1.5857027462501208e+03
-best_mse["temperature_mean"] = 6.3088676948654218e-04
-best_mse["ql_mean"] = 9.4722581733362244e+02
-best_mse["thetal_mean"] = 4.1952296308447006e-04
-best_mse["Hvar_mean"] = 8.1275664794363738e+03
-best_mse["QTvar_mean"] = 2.7043370821860535e+03
+best_mse["qt_mean"] = 1.7995063250630168e+00
+best_mse["updraft_area"] = 3.1529411391564761e+04
+best_mse["updraft_w"] = 1.0325943376733990e+03
+best_mse["updraft_qt"] = 3.1156138180627952e+01
+best_mse["updraft_thetal"] = 1.1001893832675252e+02
+best_mse["v_mean"] = 2.9275592916609139e+02
+best_mse["u_mean"] = 1.6873153880206969e+03
+best_mse["tke_mean"] = 1.3627971247685127e+03
+best_mse["temperature_mean"] = 6.8597416865395539e-04
+best_mse["ql_mean"] = 1.0669889889408917e+03
+best_mse["thetal_mean"] = 8.1757856795585861e-03
+best_mse["Hvar_mean"] = 6.8922213953686887e+03
+best_mse["QTvar_mean"] = 2.5975685843697461e+03
 
 @testset "TRMM_LBA" begin
     println("Running TRMM_LBA...")

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_

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