From 3e942184d6b62a5856f7c860d06ec5135d1d922e Mon Sep 17 00:00:00 2001
From: Anna Jaruga <ajaruga@caltech.edu>
Date: Mon, 28 Aug 2023 16:05:55 -0700
Subject: [PATCH] Dont assume that MSLP and theta ref pressure are the same

---
 src/Parameters.jl |  1 +
 src/isentropic.jl |  8 ++++----
 src/relations.jl  | 10 +++++-----
 test/relations.jl | 20 +++++++++++---------
 4 files changed, 21 insertions(+), 18 deletions(-)

diff --git a/src/Parameters.jl b/src/Parameters.jl
index 40d2f1f3..b063691d 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -21,6 +21,7 @@ param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...);
 Base.@kwdef struct ThermodynamicsParameters{FT}
     T_0::FT
     MSLP::FT
+    p_ref_theta::FT
     cp_v::FT
     cp_l::FT
     cp_i::FT
diff --git a/src/isentropic.jl b/src/isentropic.jl
index f49d57a0..7543f951 100644
--- a/src/isentropic.jl
+++ b/src/isentropic.jl
@@ -28,10 +28,10 @@ function air_pressure_given_θ(
     Φ::FT,
     ::DryAdiabaticProcess,
 ) where {FT <: AbstractFloat}
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     _R_d::FT = TP.R_d(param_set)
     _cp_d::FT = TP.cp_d(param_set)
-    return MSLP * (1 - Φ / (θ * _cp_d))^(_cp_d / _R_d)
+    return p0 * (1 - Φ / (θ * _cp_d))^(_cp_d / _R_d)
 end
 
 """
@@ -72,6 +72,6 @@ function air_temperature(
 ) where {FT <: AbstractFloat}
     _R_d::FT = TP.R_d(param_set)
     _cp_d::FT = TP.cp_d(param_set)
-    MSLP::FT = TP.MSLP(param_set)
-    return (p / MSLP)^(_R_d / _cp_d) * θ
+    p0::FT = TP.p_ref_theta(param_set)
+    return (p / p0)^(_R_d / _cp_d) * θ
 end
diff --git a/src/relations.jl b/src/relations.jl
index 38e0583d..9c6e76f3 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -2293,12 +2293,12 @@ function air_temperature_given_ρθq(
     q::PhasePartition{FT} = q_pt_0(FT),
 ) where {FT <: Real}
 
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     cvm = cv_m(param_set, q)
     cpm = cp_m(param_set, q)
     R_m = gas_constant_air(param_set, q)
     κ = 1 - cvm / cpm
-    T_u = (ρ * R_m * θ_liq_ice / MSLP)^(R_m / cvm) * θ_liq_ice
+    T_u = (ρ * R_m * θ_liq_ice / p0)^(R_m / cvm) * θ_liq_ice
     T_1 = latent_heat_liq_ice(param_set, q) / cvm
     T_2 = -κ / (2 * T_u) * (latent_heat_liq_ice(param_set, q) / cvm)^2
     return T_u + T_1 + T_2
@@ -2545,13 +2545,13 @@ function exner_given_pressure(
     p::FT,
     q::PhasePartition{FT} = q_pt_0(FT),
 ) where {FT <: Real}
-    MSLP::FT = TP.MSLP(param_set)
+    p0::FT = TP.p_ref_theta(param_set)
     # gas constant and isobaric specific heat of moist air
     _R_m = gas_constant_air(param_set, q)
     _cp_m = cp_m(param_set, q)
 
-    # return (p / MSLP)^(_R_m / _cp_m)
-    return pow_hack(p / MSLP, _R_m / _cp_m)
+    # return (p / p0)^(_R_m / _cp_m)
+    return pow_hack(p / p0, _R_m / _cp_m)
 end
 
 """
diff --git a/test/relations.jl b/test/relations.jl
index 0680ab1b..7fe55a21 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -86,7 +86,7 @@ compare_moisture(param_set, ts::PhaseNonEquil, q_pt::PhasePartition) = all((
         _T_triple = FT(TP.T_triple(param_set))
         _T_freeze = FT(TP.T_freeze(param_set))
         _T_min = FT(TP.T_min(param_set))
-        _MSLP = FT(TP.MSLP(param_set))
+        _p_ref_theta = FT(TP.p_ref_theta(param_set))
         _T_max = FT(TP.T_max(param_set))
         _kappa_d = FT(TP.kappa_d(param_set))
 
@@ -105,13 +105,13 @@ compare_moisture(param_set, ts::PhaseNonEquil, q_pt::PhasePartition) = all((
         # TODO: Use reasonable values for ambient temperature/pressure
         T∞, p∞ = T .* perturbation, p .* perturbation
         @test air_temperature.(param_set, p, θ_liq_ice, DryAdiabaticProcess()) ≈
-              (p ./ _MSLP) .^ (_R_d / _cp_d) .* θ_liq_ice
+              (p ./ _p_ref_theta) .^ (_R_d / _cp_d) .* θ_liq_ice
         @test TD.air_pressure_given_θ.(
             param_set,
             θ_liq_ice,
             Φ,
             DryAdiabaticProcess(),
-        ) ≈ _MSLP .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d)
+        ) ≈ _p_ref_theta .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d)
         @test air_pressure.(param_set, T, T∞, p∞, DryAdiabaticProcess()) ≈
               p∞ .* (T ./ T∞) .^ (FT(1) / _kappa_d)
     end
@@ -142,7 +142,7 @@ end
     _T_triple = FT(TP.T_triple(param_set))
     _T_freeze = FT(TP.T_freeze(param_set))
     _T_min = FT(TP.T_min(param_set))
-    _MSLP = FT(TP.MSLP(param_set))
+    _p_ref_theta = FT(TP.p_ref_theta(param_set))
     _T_max = FT(TP.T_max(param_set))
     _kappa_d = FT(TP.kappa_d(param_set))
     _T_icenuc = FT(TP.T_icenuc(param_set))
@@ -446,13 +446,16 @@ end
 
     # potential temperatures
     T = FT(300)
-    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _MSLP) === T
-    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _MSLP / 10) ≈
-          T * 10^(_R_d / _cp_d)
+    @test TD.liquid_ice_pottemp_given_pressure(param_set, T, _p_ref_theta) === T
     @test TD.liquid_ice_pottemp_given_pressure(
         param_set,
         T,
-        _MSLP / 10,
+        _p_ref_theta / 10,
+    ) ≈ T * 10^(_R_d / _cp_d)
+    @test TD.liquid_ice_pottemp_given_pressure(
+        param_set,
+        T,
+        _p_ref_theta / 10,
         PhasePartition(FT(1)),
     ) ≈ T * 10^(_R_v / _cp_v)
 
@@ -992,7 +995,6 @@ end
     for ArrayType in array_types
         FT = eltype(ArrayType)
         param_set = parameter_set(FT)
-        _MSLP = FT(TP.MSLP(param_set))
 
         profiles = TestedProfiles.PhaseDryProfiles(param_set, ArrayType)
         @unpack T, p, e_int, h, ρ, θ_liq_ice, phase_type = profiles