Update dry air internal energy and enthalpy definitions (eq 8a and 12a)

src/relations.jl

@@ -382,11 +382,14 @@ and, optionally,
     cvm = cv_m(param_set, q),
 ) where {FT <: Real}
     T_0 = TP.T_0(param_set)
+    R_d = TP.R_d(param_set)
     e_int_v0 = TP.e_int_v0(param_set)
     e_int_i0 = TP.e_int_i0(param_set)
     return T_0 +
            (
-        e_int - (q.tot - q.liq) * e_int_v0 + q.ice * (e_int_v0 + e_int_i0)
+        e_int - (q.tot - q.liq - q.ice) * e_int_v0 +
+        q.ice * e_int_i0 +
+        (1 - q.tot) * R_d * T_0
     ) / cvm
 end
 
@@ -407,15 +410,9 @@ and, optionally,
 ) where {FT <: Real}
     cp_m_ = cp_m(param_set, q)
     T_0 = TP.T_0(param_set)
-    R_d = TP.R_d(param_set)
     LH_v0 = TP.LH_v0(param_set)
     LH_f0 = TP.LH_f0(param_set)
-    e_int_i0 = TP.e_int_i0(param_set)
-    q_vap = vapor_specific_humidity(q)
-    return (
-        h + cp_m_ * T_0 - q_vap * LH_v0 + q.ice * LH_f0 -
-        (1 - q.tot) * R_d * T_0
-    ) / cp_m_
+    return T_0 + (h - (q.tot - q.liq - q.ice) * LH_v0 + q.ice * LH_f0) / cp_m_
 end
 
 """
@@ -469,10 +466,11 @@ and, optionally,
     cvm = cv_m(param_set, q),
 ) where {FT <: Real}
     T_0 = TP.T_0(param_set)
+    R_d = TP.R_d(param_set)
     e_int_v0 = TP.e_int_v0(param_set)
     e_int_i0 = TP.e_int_i0(param_set)
-    return cvm * (T - T_0) + (q.tot - q.liq) * e_int_v0 -
-           q.ice * (e_int_v0 + e_int_i0)
+    return cvm * (T - T_0) + (q.tot - q.liq - q.ice) * e_int_v0 -
+           q.ice * e_int_i0 - (1 - q.tot) * R_d * T_0
 end
 
 """
@@ -516,8 +514,8 @@ The dry air internal energy
 @inline function internal_energy_dry(param_set::APS, T::FT) where {FT <: Real}
     T_0 = TP.T_0(param_set)
     cv_d = TP.cv_d(param_set)
-
-    return cv_d * (T - T_0)
+    R_d = TP.R_d(param_set)
+    return cv_d * (T - T_0) - R_d * T_0
 end
 
 """
@@ -1331,7 +1329,7 @@ is a function that is 1 above `T_freeze` and goes to zero below `T_icenuc`.
     # Model for cloud water condensate probabilty when temperature is below
     # freezing (all liquid condensate), but above the temperature of homogeneous
     # ice nucleation (all ice condensate). In it's simplest form, this is simply
-    # a number that varies between 0 and 1. For example, see figure 6 in 
+    # a number that varies between 0 and 1. For example, see figure 6 in
     # Hu et al., https://doi.org/10.1029/2009JD012384, JGR (2010).
     λᵖ = ((T - Tⁿ) / (Tᶠ - Tⁿ))^α
 

test/relations.jl

@@ -221,7 +221,7 @@ end
     # test the wrapper for q_vap_saturation over liquid water and ice
     ρ = FT(1)
     ρu = FT[1, 2, 3]
-    ρe = FT(1100)
+    ρe = FT(1100) - _R_d * _T_triple
     e_pot = FT(93)
     e_int = internal_energy(ρ, ρe, ρu, e_pot)
     q_pt = PhasePartition(FT(0.02), FT(0.002), FT(0.002))
@@ -291,31 +291,32 @@ end
     T = FT(300)
     e_kin = FT(11)
     e_pot = FT(13)
-    @test air_temperature(param_set, _cv_d * (T - _T_0)) === FT(T)
+    @test air_temperature(param_set, _cv_d * (T - _T_0) - _R_d * _T_triple) ===
+          FT(T)
     @test air_temperature(
         param_set,
-        _cv_d * (T - _T_0),
+        _cv_d * (T - _T_0) - _R_d * _T_triple,
         PhasePartition(FT(0)),
     ) === FT(T)
 
     @test air_temperature(
         param_set,
-        cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0),
+        cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0) - _R_d * _T_triple,
         PhasePartition(FT(0)),
     ) === FT(T)
     @test air_temperature(
         param_set,
-        cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) +
-        q_tot * _e_int_v0,
+        cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) -
+        (1 - q_tot) * _R_d * _T_triple + q_tot * _e_int_v0,
         PhasePartition(q_tot),
     ) ≈ FT(T)
 
     @test total_energy(param_set, FT(e_kin), FT(e_pot), _T_0) ===
-          FT(e_kin) + FT(e_pot)
+          FT(e_kin) + FT(e_pot) - _R_d * _T_triple
     @test total_energy(param_set, FT(e_kin), FT(e_pot), FT(T)) ≈
-          FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0)
+          FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0) - _R_d * _T_triple
     @test total_energy(param_set, FT(0), FT(0), _T_0, PhasePartition(q_tot)) ≈
-          q_tot * _e_int_v0
+          q_tot * _e_int_v0 - (1 - q_tot) * _R_d * _T_triple
 
     # phase partitioning in equilibrium
     q_liq = FT(0.1)
@@ -595,7 +596,13 @@ end
             )
         _e_int = (e_int_upper .+ e_int_lower) / 2
         ts = PhaseEquil_ρeq.(param_set, ρ, _e_int, q_tot)
-        @test all(air_temperature.(param_set, ts) .== Ref(_T_freeze))
+        @test all(
+            isapprox.(
+                air_temperature.(param_set, ts),
+                Ref(_T_freeze),
+                rtol = rtol_temperature,
+            ),
+        )
 
         # Args needs to be in sync with PhaseEquil:
         ts =
@@ -608,7 +615,13 @@ end
                 FT(1e-1),
                 RS.SecantMethod,
             )
-        @test all(air_temperature.(param_set, ts) .== Ref(_T_freeze))
+        @test all(
+            isapprox.(
+                air_temperature.(param_set, ts),
+                Ref(_T_freeze),
+                rtol = rtol_temperature,
+            ),
+        )
 
         # PhaseEquil
         ts_exact = PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot, 100, FT(1e-6))