Add quadratic term to air temp expansion

src/Common/MoistThermodynamics/MoistThermodynamics.jl

@@ -927,9 +927,14 @@ and, optionally,
 function air_temperature_from_liquid_ice_pottemp(θ_liq_ice::FT, ρ::FT,
   q::PhasePartition{FT}=q_pt_0(FT)) where {FT<:Real}
 
-  cv = cv_m(q)
+  cvm = cv_m(q)
+  cpm = cp_m(q)
   R_m = gas_constant_air(q)
-  return θ_liq_ice * (ρ*R_m*θ_liq_ice/FT(MSLP))^(R_m/cv) + latent_heat_liq_ice(q)/cv
+  κ = 1 - cvm/cpm
+  T_u = (ρ*R_m*θ_liq_ice/FT(MSLP))^(R_m/cvm)*θ_liq_ice
+  T_1 = latent_heat_liq_ice(q)/cvm
+  T_2 = -κ/(2*T_u)*(latent_heat_liq_ice(q)/cvm)^2
+  return T_u + T_1 + T_2
 end
 
 """

test/Common/MoistThermodynamics/runtests.jl

@@ -237,6 +237,8 @@ end
 
     # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice)
     T_non_linear = air_temperature_from_liquid_ice_pottemp_non_linear.(θ_liq_ice, ρ, FT(1e-3), 10, q_pt)
+    T_expansion = air_temperature_from_liquid_ice_pottemp.(θ_liq_ice, ρ, q_pt)
+    @test all(isapprox.(T_non_linear, T_expansion, rtol=rtol))
     e_int_ = internal_energy.(T_non_linear, q_pt)
     ts = PhaseNonEquil.(e_int_, ρ, q_pt)
     @test all(T_non_linear .≈ air_temperature.(ts))