From d8ab91a320202ad9f908ff1851036926e226bb2e Mon Sep 17 00:00:00 2001
From: Anna Jaruga <ajaruga@caltech.edu>
Date: Thu, 3 Mar 2022 14:51:20 -0800
Subject: [PATCH] Add more flexible phase partition in equilibrium

---
 Project.toml        |  4 ++--
 gpuenv/Project.toml |  2 +-
 src/relations.jl    | 28 +++++++++++++++++++++-------
 test/Project.toml   |  2 +-
 test/relations.jl   |  2 +-
 5 files changed, 26 insertions(+), 12 deletions(-)

diff --git a/Project.toml b/Project.toml
index 80d0925c..221c6f64 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Thermodynamics"
 uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 authors = ["Climate Modeling Alliance"]
-version = "0.5.11"
+version = "0.6.0"
 
 [deps]
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
@@ -11,7 +11,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 
 [compat]
-CLIMAParameters = "0.1, 0.2, 0.3, 0.4"
+CLIMAParameters = "0.4"
 DocStringExtensions = "0.8.1"
 KernelAbstractions = "0.7.2"
 RootSolvers = "0.2, 0.3"
diff --git a/gpuenv/Project.toml b/gpuenv/Project.toml
index 6c53745a..aa62f11d 100644
--- a/gpuenv/Project.toml
+++ b/gpuenv/Project.toml
@@ -10,7 +10,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
-CLIMAParameters = "0.1, 0.2, 0.3, 0.4"
+CLIMAParameters = "0.4"
 CUDA = "3.5"
 CUDAKernels = "0.3"
 KernelAbstractions = "0.7.2"
diff --git a/src/relations.jl b/src/relations.jl
index 5d21dc88..0ab06702 100644
--- a/src/relations.jl
+++ b/src/relations.jl
@@ -1163,7 +1163,7 @@ them.
 
 # ThermodynamicState
 Otherwise, phase equilibrium is assumed so that the fraction of liquid
-is a function that is 1 above `T_freeze` and goes to zero below `T_freeze`.
+is a function that is 1 above `T_freeze` and goes to zero below `T_icenuc`.
 """
 function liquid_fraction(
     param_set::APS,
@@ -1172,7 +1172,16 @@ function liquid_fraction(
     q::PhasePartition{FT} = q_pt_0(FT),
 ) where {FT <: Real, phase_type <: ThermodynamicState}
     _T_freeze::FT = CPP.T_freeze(param_set)
-    return FT(T > _T_freeze)
+    _T_icenuc::FT = CPP.T_icenuc(param_set)
+    _pow_icenuc::FT = CPP.pow_icenuc(param_set)
+
+    if T > _T_freeze
+        return FT(1)
+    elseif (T > _T_icenuc && T <= _T_freeze)
+        return ((T - _T_icenuc) / (_T_freeze - _T_icenuc))^_pow_icenuc
+    else
+        return FT(0)
+    end
 end
 
 function liquid_fraction(
@@ -1293,19 +1302,24 @@ function ∂e_int_∂T(
     cv_i::FT = CPP.cv_i(param_set)
     e_int_v0::FT = CPP.e_int_v0(param_set)
     e_int_i0::FT = CPP.e_int_i0(param_set)
+    T_f::FT = CPP.T_freeze(param_set)
+    T_i::FT = CPP.T_icenuc(param_set)
+    n_i::FT = CPP.pow_icenuc(param_set)
 
-    cvm = cv_m(
-        param_set,
-        PhasePartition_equil(param_set, T, ρ, q_tot, phase_type),
-    )
+    q = PhasePartition_equil(param_set, T, ρ, q_tot, phase_type)
+    q_c = condensate(q)
+    cvm = cv_m(param_set, q)
     q_vap_sat = q_vap_saturation(param_set, T, ρ, phase_type)
     λ = liquid_fraction(param_set, T, phase_type)
     L = λ * LH_v0 + (1 - λ) * LH_s0
+
+    ∂λ_∂T = (T_i < T < T_f) ? (1 / (T_f - T_i))^n_i * n_i * T^(n_i - 1) : FT(0)
     ∂q_vap_sat_∂T = q_vap_sat * L / (R_v * T^2)
     dcvm_dq_vap = cv_v - λ * cv_l - (1 - λ) * cv_i
     return cvm +
            (e_int_v0 + (1 - λ) * e_int_i0 + (T - T_0) * dcvm_dq_vap) *
-           ∂q_vap_sat_∂T
+           ∂q_vap_sat_∂T +
+           q_c * e_int_i0 * ∂λ_∂T
 end
 
 """
diff --git a/test/Project.toml b/test/Project.toml
index eb95a521..b8d8c768 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -20,7 +20,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
-CLIMAParameters = "0.1, 0.2, 0.3, 0.4"
+CLIMAParameters = "0.4"
 KernelAbstractions = "0.7.2"
 RootSolvers = "0.2, 0.3"
 UnPack = "1.0"
diff --git a/test/relations.jl b/test/relations.jl
index 9194ef16..97889700 100644
--- a/test/relations.jl
+++ b/test/relations.jl
@@ -703,7 +703,7 @@ end
 
         @test count(air_temperature.(ts_lower) .== Ref(_T_freeze)) ≥ 217
         @test count(air_temperature.(ts_upper) .== Ref(_T_freeze)) ≥ 217
-        @test count(air_temperature.(ts_mid) .== Ref(_T_freeze)) ≥ 1395
+        @test count(air_temperature.(ts_mid) .== Ref(_T_freeze)) ≥ 1392
         # we should do this instead, but we're failing because some inputs are bad
         # E.g. p ~ 110_000 Pa, q_tot ~ 0.16, which results in negative θ_liq_ice
         # This means that we should probably update our tested profiles.