From 60f82885ff5da339032155b491f6e65ffc3ca1fa Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 12 Sep 2024 15:36:25 -0700 Subject: [PATCH] Update dry air internal energy and enthalpy definitions (eq 8a and 12a) --- src/relations.jl | 24 +++++++++++------------- test/relations.jl | 35 ++++++++++++++++++++++++----------- 2 files changed, 35 insertions(+), 24 deletions(-) diff --git a/src/relations.jl b/src/relations.jl index 721bf11b..e69ebef6 100644 --- a/src/relations.jl +++ b/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ⁿ))^α diff --git a/test/relations.jl b/test/relations.jl index 9bda6181..177fdbd4 100644 --- a/test/relations.jl +++ b/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))