From 7a578b6571b92e6c6a8dbb33aaf97c2278131431 Mon Sep 17 00:00:00 2001
From: Haakon Ludvig Langeland Ervik
 <45243236+haakon-e@users.noreply.github.com>
Date: Mon, 22 Jun 2026 16:53:27 -0700
Subject: [PATCH] fix: units of Barklie Gokhale 1959

See: https://github.com/CliMA/ClimaParams.jl/pull/295
---
 Project.toml                    |  2 +-
 src/IceNucleation.jl            |  7 +++----
 src/parameters/IceNucleation.jl | 18 ++++++++----------
 test/gpu_tests.jl               |  2 +-
 4 files changed, 13 insertions(+), 16 deletions(-)

diff --git a/Project.toml b/Project.toml
index 6cd35a786..65b86c203 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7"
 EmulatorModelsExt = ["DataFrames", "MLJ"]
 
 [compat]
-ClimaParams = "1"
+ClimaParams = "1.0.18"
 DataFrames = "1.6"
 DocStringExtensions = "0.8, 0.9"
 FastGaussQuadrature = "1"
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
index bda69cda8..ff4bdbc36 100644
--- a/src/IceNucleation.jl
+++ b/src/IceNucleation.jl
@@ -172,8 +172,8 @@ Compute number of ice crystals formed from heterogeneous condensation freezing
 
 # Arguments
  - `ip`: The [`CMP.MorrisonMilbrandt2014`](@ref) paramterization, where:
-    + `het_a`: empirical parameter [C⁻¹]
-    + `het_B`: water-type dependent parameter [cm⁻³ s⁻¹]
+    + `het_a`: empirical parameter [K⁻¹]
+    + `het_B`: water-type dependent parameter [m⁻³ s⁻¹]
     + `T₀`: freezing temperature [K]
  - `T`: air temperature [K],
  - `Nₗ`: number of droplets [m⁻³],
@@ -200,9 +200,8 @@ Nᵤ(t) = N₀ exp(-B Vₗ [exp(aTₛ)] t)
 ```
 """
 function P3_het_N_i((; het_a, het_B, T₀)::CMP.MorrisonMilbrandt2014, T, Nₗ, Vₗ, Δt)
-    Vₗ_cm³ = Vₗ * 1_000_000  # converted from m^3 to cm^3
     Tₛ = T₀ - T
-    return Nₗ * (1 - exp(-het_B * Vₗ_cm³ * Δt * exp(het_a * Tₛ)))
+    return Nₗ * (1 - exp(-het_B * Vₗ * Δt * exp(het_a * Tₛ)))
 end
 
 """
diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl
index ca070dc15..24c7a38ad 100644
--- a/src/parameters/IceNucleation.jl
+++ b/src/parameters/IceNucleation.jl
@@ -87,9 +87,9 @@ $(DocStringExtensions.FIELDS)
     c₂::FT
     "T₀"
     T₀::FT
-    "heterogeneous freezing parameter a [°C^-1]"
+    "heterogeneous freezing parameter a [K⁻¹]"
     het_a::FT
-    "heterogeneous freezing parameter B [cm^-3 s^-1]"
+    "heterogeneous freezing parameter B [m⁻³ s⁻¹]"
     het_B::FT
 end
 
@@ -121,30 +121,28 @@ $(DocStringExtensions.FIELDS)
 
 # Callable interface
 
-    (rf::RainFreezing)(T, T₀) → het_B * 10⁶ * exp(het_a * (T₀ - T))
+    (rf::RainFreezing)(T, T₀) → het_B * exp(het_a * (T₀ - T))
 
-Returns the volumetric freezing rate in SI units [m⁻³ s⁻¹].
-The stored `het_B` is in [cm⁻³ s⁻¹]; the factor 10⁶ converts cm³ → m³.
+Compute the volumetric freezing rate [m⁻³ s⁻¹]
 """
 @kwdef struct RainFreezing{FT} <: ParametersType
-    "empirical parameter [°K⁻¹]"
+    "empirical parameter [K⁻¹]"
     het_a::FT
-    "water-type dependent parameter [cm⁻³ s⁻¹]"
+    "water-type dependent parameter [m⁻³ s⁻¹]"
     het_B::FT
 end
 
 function RainFreezing(td::CP.ParamDict)
     name_map = (;
         :BarklieGokhale1959_a_parameter => :het_a,
-        :BarklieGokhale1959_B_parameter => :het_B,  # TODO: Fix units in CP, then here
+        :BarklieGokhale1959_B_parameter => :het_B,
     )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     return RainFreezing(; parameters...)
 end
 
 # Callable: returns the Bigg (1953) volumetric freezing rate [m⁻³(water) s⁻¹]
-# het_B is stored in [cm⁻³ s⁻¹]; multiply by 10⁶ to convert to [m⁻³(water) s⁻¹].
-(rf::RainFreezing)(T, T₀) = rf.het_B * 1_000_000 * exp(rf.het_a * (T₀ - T))
+(rf::RainFreezing)(T, T₀) = rf.het_B * exp(rf.het_a * (T₀ - T))
 
 """
     IceNucleationParameters{FT, DEP, HOM, P3_type}
diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl
index 7477357ca..f5b54e2db 100644
--- a/test/gpu_tests.jl
+++ b/test/gpu_tests.jl
@@ -1003,7 +1003,7 @@ function test_gpu(FT)
                 FT(0.0002736160475969029)
             else
                 # loss of precision due to
-                # `exp(-B * V_l_converted * Δt * exp(a * Tₛ))` -> 0.9999999 (Float32)
+                # `exp(-het_B * V_l * Δt * exp(het_a * Tₛ))` -> 0.9999999 (Float32)
                 # instead of 0.9999998631919762 (Float64).
                 FT(0.00023841858f0)
             end