From 4e03889f606e2750371280a8382c4743553a926b Mon Sep 17 00:00:00 2001
From: Anastasia Popova <apopova@caltech.edu>
Date: Thu, 1 Feb 2024 13:36:49 -0800
Subject: [PATCH] unitless N_hat added, Float32 tests passing

---
 docs/src/P3Scheme.md | 2 +-
 src/P3Scheme.jl      | 6 ++++--
 test/p3_tests.jl     | 8 ++++----
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md
index 83f9176e9..0c8ffbbe2 100644
--- a/docs/src/P3Scheme.md
+++ b/docs/src/P3Scheme.md
@@ -91,7 +91,7 @@ N'(D) = N_{0} D^\mu \, e^{-\lambda \, D}
 where:
  - ``N'`` is the number concentration in ``m^{-4}``
  - ``D`` is the maximum particle dimension in ``m``,
- - ``N_0`` is the intercept parameter in ``m^{-4}``,
+ - ``N_0`` is the intercept parameter in ``m^{-5 - \mu }``,
  - ``\mu`` is the shape parameter (dimensionless),
  - ``\lambda`` is the slope parameter in ``m^{-1}``.
 
diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl
index f89ab8c50..15b1eb5c0 100644
--- a/src/P3Scheme.jl
+++ b/src/P3Scheme.jl
@@ -125,7 +125,7 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT}
     @assert F_r < FT(1)    # ... and there must always be some unrimed part
 
     if F_r == FT(0)
-        return (; D_cr = 0, D_gr = 0, ρ_g = 0, ρ_d = 0)
+        return (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0))
     else
         @assert ρ_r > FT(0)   # rime density must be positive ...
         @assert ρ_r <= p3.ρ_l # ... and as a bulk ice density can't exceed the density of water
@@ -344,7 +344,9 @@ function distribution_parameter_solver(
     th = thresholds(p3, ρ_r, F_r)
 
     # To ensure that λ is positive solve for x such that λ = exp(x)
-    shape_problem(x) = q - q_gamma(p3, F_r, N, x, th)
+    # We divide by N̂ to deal with large N₀ values for Float32
+    N̂ = FT(1e30)
+    shape_problem(x) = q / N̂ - q_gamma(p3, F_r, N / N̂, x, th)
 
     # Find slope parameter
     x =
diff --git a/test/p3_tests.jl b/test/p3_tests.jl
index 052e2f95f..c5a0ecd68 100644
--- a/test/p3_tests.jl
+++ b/test/p3_tests.jl
@@ -156,7 +156,7 @@ function test_p3_shape_solver(FT)
     TT.@testset "shape parameters (nonlinear solver function)" begin
 
         # initialize test values:
-        eps = FT(1e-3)
+        ep = 1e3 * eps(FT)
         N_test = (FT(1e8))                             # N values
         λ_test = (FT(15000), FT(20000))                # test λ values in range
         ρ_r_test = (FT(200)) #, FT(1)) #, FT(100))     # representative ρ_r values
@@ -188,8 +188,8 @@ function test_p3_shape_solver(FT)
                         )
 
                         # Compare solved values with the input expected values
-                        TT.@test λ ≈ λ_ex rtol = eps
-                        TT.@test N₀ ≈ N₀_ex rtol = eps
+                        TT.@test λ ≈ λ_ex rtol = ep
+                        TT.@test N₀ ≈ N₀_ex rtol = ep
                     end
                 end
             end
@@ -202,7 +202,7 @@ test_p3_thresholds(Float32)
 test_p3_mass(Float32)
 #TODO - only works for Float64 now. We should switch the units inside the solver
 # from SI base to something more managable
-#test_p3_shape_solver(Float32)
+test_p3_shape_solver(Float32)
 
 println("Testing Float64")
 test_p3_thresholds(Float64)