From 08704d24cf3d9ef10bcf3fb04a97a65cba8385dc Mon Sep 17 00:00:00 2001
From: amylu00 <alu3@caltech.edu>
Date: Mon, 8 Apr 2024 17:52:39 -0700
Subject: [PATCH] derived \Delta a_w equivalent parameterization

---
 docs/src/plots/linear_HOM_J.jl | 52 ++++++++++++++++++----------------
 parcel/ParcelTendencies.jl     |  2 +-
 src/IceNucleation.jl           |  4 +--
 3 files changed, 31 insertions(+), 27 deletions(-)

diff --git a/docs/src/plots/linear_HOM_J.jl b/docs/src/plots/linear_HOM_J.jl
index 9f3b3349f5..3e425d2d1e 100644
--- a/docs/src/plots/linear_HOM_J.jl
+++ b/docs/src/plots/linear_HOM_J.jl
@@ -10,52 +10,56 @@ const CMO = CM.Common
 FT = Float32
 const tps = TD.Parameters.ThermodynamicsParameters(FT)
 
-# From Atkinson, Murray, and O'Sullivan 2016 (DOI:10.1021/acs.jpca.6b03843)
+# From Atkinson, Murray, and O'Sullivan 2016
+# (DOI:10.1021/acs.jpca.6b03843)
+# Using overall volume fit
 function Atkinson_J(T)
     lnJ = -4.0106 * T + 963.706
-    return exp(lnJ) * 1e6
+    return exp(lnJ)
 end
 
 # Initializing
-T_range = range(FT(229.0), stop = 238, length = 50)  # air temperature
+T_range = range(FT(234.0), stop = 237, length = 50)  # air temperature
 Sₗ = FT(1.1)
 eₛ = [TD.saturation_vapor_pressure(tps, T, TD.Liquid()) for T in T_range]
-e = Sₗ .* eₛ
+e_range = Sₗ .* eₛ
 
-Δa = zeros(FT, 1, length(T_range))
+Δa = zeros(FT, length(T_range))
 
 # Solving for Δa and J values
-J = [Atkinson_J(T) for T in T_range]
+J = [Atkinson_J(T) for T in T_range] # [cm^-3 s^-1]
 log10J = @. log10(J)
-Δa = [CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T) for T in T_range]
+for (i, T) in enumerate(T_range)
+    Δa[i] = CMO.a_w_eT(tps, e_range[i], T) - CMO.a_w_ice(tps, T)
+end
+
+# Obtaining a linear parameterization for J in terms of Δa
+slope = (log10J[end] - log10J[1]) / (Δa[end] - Δa[1])
+intercept = slope * (Δa[1]) + log10J[1]
 
 # Plotting J vs Δa
-fig = MK.Figure(resolution = (800, 500))
+fig = MK.Figure(resolution = (800, 700))
 ax1 = MK.Axis(
     fig[1, 1],
-    ylabel = "log10(J) with J in SI units",
+    ylabel = "log10(J), J = [cm^-3 s^-1]",
     xlabel = "Temperature [K]",
     title = "log10(J) vs T",
-    xticklabelsize = 14.0f0,
-    xlabelsize = 14,
-    ylabelsize = 14,
-    #limits = ((228.0, 240.0), nothing),
 )
 ax2 = MK.Axis(
-    fig[1, 2],
+    fig[2, 1],
     xlabel = "Δa_w [-]",
-    ylabel = "log10(J) [cm^-3 s^-1]",
+    ylabel = "log10(J), J = [cm^-3 s^-1]",
     title = "log10(J) vs Δa_w",
-    xticklabelsize = 14.0f0,
-    xlabelsize = 14,
-    ylabelsize = 14,
 )
 
-MK.lines!(ax1, T_range, J)
-MK.lines!(ax2, Δa, J)
-#! format: on
+MK.lines!(ax1, T_range, log10J, label = "ln(J) = -4.0106 T + 963.706",)
+MK.lines!(ax2, Δa, log10J,
+    label = "log10(J) = $slope Δa + $intercept",
+)
+
+MK.axislegend(ax1, position = :rt, labelsize = 18.0f0)
+MK.axislegend(ax2, position = :rb, labelsize = 18.0f0)
 
-MK.axislegend(ax1, position = :rt, labelsize = 13.0f0)
-MK.axislegend(ax2, position = :rb, labelsize = 14.0f0)
+#! format: on
 
-MK.save("HomLinearJ.svg", fig)
+MK.save("linear_HOM_J.svg", fig)
diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl
index 84d4630abc..245854be5e 100644
--- a/parcel/ParcelTendencies.jl
+++ b/parcel/ParcelTendencies.jl
@@ -180,7 +180,7 @@ function homogeneous_freezing(params::ABHOM, PSD, state)
         Δa_w = ips.homogeneous.Δa_w_max
     end
 
-    J = CMI_hom.homogeneous_J(ips.homogeneous, Δa_w)
+    J = CMI_hom.homogeneous_J_cubic(ips.homogeneous, Δa_w)
 
     return min(max(FT(0), J * Nₗ * PSD.Vₗ), Nₗ / const_dt)
 end
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
index 51be1ccbc2..f2ad1a1924 100644
--- a/src/IceNucleation.jl
+++ b/src/IceNucleation.jl
@@ -247,9 +247,9 @@ Returns the homogeneous freezing nucleation rate coefficient,
 `J`, in m^-3 s^-1 for sulphuric acid solutions.
 Parameterization based on Atkinson et al 2016, DOI:10.1021/acs.jpca.6b03847.
 """
-function homogeneous_J_linear(ip::CMP.Koop2000, Δa_w::FT) where {FT}
+function homogeneous_J_linear(ip::CMP.Atkinson2016, Δa_w::FT) where {FT}
 
-    logJ::FT = ip.c₁ + ip.c₂ * Δa_w
+    logJ::FT = ip.c₁ * Δa_w + ip.c₂
 
     return FT(10)^(logJ) * 1e6
 end