diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 3a4a6600a..d5ca82570 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1,3 +1,5 @@
+margin = 80
+indent = 4
 always_for_in = true
 whitespace_typedefs = true
-whitespace_ops_in_indices = true
+whitespace_ops_in_indices = true
\ No newline at end of file
diff --git a/.dev/up_deps.jl b/.dev/up_deps.jl
index cb60f141f..69190000f 100644
--- a/.dev/up_deps.jl
+++ b/.dev/up_deps.jl
@@ -4,7 +4,12 @@ files in all of our environments.
 =#
 
 root = dirname(@__DIR__)
-dirs = (root, joinpath(root, "test"), joinpath(root, "docs"), joinpath(root, "parcel"))
+dirs = (
+    root,
+    joinpath(root, "test"),
+    joinpath(root, "docs"),
+    joinpath(root, "parcel"),
+)
 
 cd(root) do
     for dir in dirs
diff --git a/box/Alpert_Knopf_2016_forward.jl b/box/Alpert_Knopf_2016_forward.jl
index 97eeb3309..b489e6a2e 100644
--- a/box/Alpert_Knopf_2016_forward.jl
+++ b/box/Alpert_Knopf_2016_forward.jl
@@ -49,8 +49,12 @@ p_def = (; tps, A_aero, aerosol, cooling_rate, N₀)
 
 # run the box model without and with distribution sampling
 sol_cst = run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = false))
-sol_var =
-    run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = true, Aj_sorted = Aj_sorted))
+sol_var = run_box(
+    IC,
+    t_0,
+    t_end,
+    (p_def..., const_dt = dt, flag = true, Aj_sorted = Aj_sorted),
+)
 
 # Compute the immersion freezing rate for the no-sampling case
 Δa = @. FT(1) - CMO.a_w_ice(tps, sol_cst[1, :])
diff --git a/docs/make.jl b/docs/make.jl
index 4b7b61e72..a00c67c6d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -72,7 +72,12 @@ pages = Any[
 ]
 
 mathengine = MathJax(
-    Dict(:TeX => Dict(:equationNumbers => Dict(:autoNumber => "AMS"), :Macros => Dict())),
+    Dict(
+        :TeX => Dict(
+            :equationNumbers => Dict(:autoNumber => "AMS"),
+            :Macros => Dict(),
+        ),
+    ),
 )
 
 format = Documenter.HTML(
diff --git a/docs/src/guides/GettingStarted.jl b/docs/src/guides/GettingStarted.jl
index dd3875373..37e9fc5e3 100644
--- a/docs/src/guides/GettingStarted.jl
+++ b/docs/src/guides/GettingStarted.jl
@@ -35,5 +35,6 @@ nothing #hide
 
 # Finally, we call `accretion`, which will return the accretion rates for
 # cloud and rain water specific humidities, as well as cloud and rain water number concentrations.
-(; dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt) = CM2.accretion(SB2006, qₗ, qᵣ, ρₐ, Nₗ)
+(; dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt) =
+    CM2.accretion(SB2006, qₗ, qᵣ, ρₐ, Nₗ)
 @info("Accretion rates: ", dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt)
diff --git a/docs/src/guides/Parameters.jl b/docs/src/guides/Parameters.jl
index ccf187b1b..3d1dcb12c 100644
--- a/docs/src/guides/Parameters.jl
+++ b/docs/src/guides/Parameters.jl
@@ -46,8 +46,10 @@ nothing #hide
 # Overwriting the parameters can also be done using dictionaries instead of toml files.
 # As an example we create a dictionary were we define a different value of the
 # rain autoconversion time scale and create another parameter struct based on it.
-override_file =
-    Dict("rain_autoconversion_timescale" => Dict("value" => 13, "type" => "float"))
+override_file = Dict(
+    "rain_autoconversion_timescale" =>
+        Dict("value" => 13, "type" => "float"),
+)
 toml_dict2 = CP.create_toml_dict(FT; override_file)
 const overwrite2 = CMP.Rain(toml_dict2)
 
diff --git a/docs/src/plots/ARGplots.jl b/docs/src/plots/ARGplots.jl
index 585eaa2e8..06d2b784c 100644
--- a/docs/src/plots/ARGplots.jl
+++ b/docs/src/plots/ARGplots.jl
@@ -48,7 +48,9 @@ function mass2vol(mass_mixing_ratios)
     else
         densities = (sulfate.ρ,)
     end
-    volfractions = (mass_mixing_ratios ./ densities) ./ sum(mass_mixing_ratios ./ densities)
+    volfractions =
+        (mass_mixing_ratios ./ densities) ./
+        sum(mass_mixing_ratios ./ densities)
     return volfractions
 end
 
@@ -160,9 +162,11 @@ function make_ARG_figX(X)
                 end
                 AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
                 act_frac1[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
+                    N_1
                 act_frac2[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N2i
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
+                    N2i
                 global it += 1
             end
 
@@ -211,9 +215,11 @@ function make_ARG_figX(X)
                 end
                 AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
                 act_frac1[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
+                    N_1
                 act_frac2[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N2i
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
+                    N2i
                 global it += 1
             end
 
@@ -264,9 +270,11 @@ function make_ARG_figX(X)
                 end
                 AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
                 act_frac1[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
+                    N_1
                 act_frac2[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N_2
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
+                    N_2
                 global it += 1
             end
 
@@ -314,9 +322,11 @@ function make_ARG_figX(X)
                 end
                 AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
                 act_frac1[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
+                    N_1
                 act_frac2[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N_2
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
+                    N_2
                 global it += 1
             end
 
@@ -366,9 +376,11 @@ function make_ARG_figX(X)
 
             for wi in w
                 act_frac1[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[1] / N_1
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[1] /
+                    N_1
                 act_frac2[it] =
-                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[2] / N_2
+                    AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[2] /
+                    N_2
                 global it += 1
             end
 
diff --git a/docs/src/plots/ARGplots_fig1.jl b/docs/src/plots/ARGplots_fig1.jl
index 89158eb16..51e709c74 100644
--- a/docs/src/plots/ARGplots_fig1.jl
+++ b/docs/src/plots/ARGplots_fig1.jl
@@ -68,7 +68,8 @@ for N_2 in N_2_range
         (sulfate.ρ,),
     )
     AD_B = AM.AerosolDistribution((paper_mode_1_B, paper_mode_2_B))
-    N_act_frac_B[it] = AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / N_1
+    N_act_frac_B[it] =
+        AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / N_1
     global it += 1
 end
 
@@ -83,7 +84,17 @@ PL.plot(
     xlabel = "Mode 2 aerosol number concentration [1/cm3]",
     ylabel = "Mode 1 number fraction activated",
 )
-PL.scatter!(Fig1_x_obs, Fig1_y_obs, markercolor = :black, label = "paper observations")
-PL.plot!(Fig1_x_param, Fig1_y_param, linecolor = :black, label = "paper parameterization")
+PL.scatter!(
+    Fig1_x_obs,
+    Fig1_y_obs,
+    markercolor = :black,
+    label = "paper observations",
+)
+PL.plot!(
+    Fig1_x_param,
+    Fig1_y_param,
+    linecolor = :black,
+    label = "paper parameterization",
+)
 
 PL.savefig("Abdul-Razzak_and_Ghan_fig_1.svg")
diff --git a/docs/src/plots/Baumgartner2022_fig5.jl b/docs/src/plots/Baumgartner2022_fig5.jl
index 7568882ef..02a59f391 100644
--- a/docs/src/plots/Baumgartner2022_fig5.jl
+++ b/docs/src/plots/Baumgartner2022_fig5.jl
@@ -35,14 +35,17 @@ sol = RS.find_zero(fun, RS.SecantMethod(FT(0), FT(0.5)), RS.CompactSolution())
 a_w = Δa_crit .+ a_w_ice
 
 # various ways Baumgartner calculates vapor pressures
-Baumgartner_p_ice(T) = exp(9.550426 - 5723.265 / T + 3.53068 * log(T) - 0.00728332 * T)
+Baumgartner_p_ice(T) =
+    exp(9.550426 - 5723.265 / T + 3.53068 * log(T) - 0.00728332 * T)
 MK_p_liq(T) = exp(
     54.842763 - 6763.22 / T - 4.21 * log(T) +
     0.000367 * T +
-    tanh(0.0415 * (T - 218.8)) * (53.878 - 1331.22 / T - 9.44523 * log(T) + 0.014025 * T),
+    tanh(0.0415 * (T - 218.8)) *
+    (53.878 - 1331.22 / T - 9.44523 * log(T) + 0.014025 * T),
+)
+Tab_p_liq(T) = exp(
+    100 * (18.4524 - 3505.15788 / T - 330918.55 / (T^2) + 12725068.26 / (T^3)),
 )
-Tab_p_liq(T) =
-    exp(100 * (18.4524 - 3505.15788 / T - 330918.55 / (T^2) + 12725068.26 / (T^3)))
 Nach_p_liq(T) = exp(74.8727 - 7167.405 / T - 7.77107 * log(T) + 0.00505 * T)
 
 a_w_MK = [Δa_crit + (Baumgartner_p_ice(T)) / (MK_p_liq(T)) for T in T_range]
@@ -50,7 +53,9 @@ a_w_Tab = [Δa_crit + (Baumgartner_p_ice(T)) / (Tab_p_liq(T)) for T in T_range]
 a_w_Nach = [Δa_crit + (Baumgartner_p_ice(T)) / (Nach_p_liq(T)) for T in T_range]
 a_w_Luo = [
     Δa_crit +
-    (Baumgartner_p_ice(T)) / (CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T)) for T in T_range
+    (Baumgartner_p_ice(T)) /
+    (CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T)) for
+    T in T_range
 ]
 
 #! format: off
diff --git a/docs/src/plots/Cholette2019_fig1.jl b/docs/src/plots/Cholette2019_fig1.jl
index 1784adcee..16cce4136 100644
--- a/docs/src/plots/Cholette2019_fig1.jl
+++ b/docs/src/plots/Cholette2019_fig1.jl
@@ -67,14 +67,20 @@ function fig1()
     F_liqs = (FT(0), FT(0.33), FT(0.67), FT(1))
     F_rs = (FT(0), FT(1 - eps(FT)))
 
-    labels = ["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1", "ϕᵢ = 1"]
+    labels =
+        ["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1", "ϕᵢ = 1"]
     colors = [:black, :blue, :red, :green]
     res = 50
 
     fig = Plt.Figure(size = (1200, 400))
 
     #F_r = 0
-    ax1 = Plt.Axis(fig[1:7, 1:9], title = "Fᵣ = 0", xlabel = "Dₚ (m)", ylabel = "V (m s⁻¹)")
+    ax1 = Plt.Axis(
+        fig[1:7, 1:9],
+        title = "Fᵣ = 0",
+        xlabel = "Dₚ (m)",
+        ylabel = "V (m s⁻¹)",
+    )
 
     for i in 1:4
         D_ps, V_0, V_0_ϕ = get_values(
@@ -87,12 +93,23 @@ function fig1()
             res,
         )
         Plt.lines!(ax1, D_ps, V_0_ϕ, color = colors[i], label = labels[i])
-        Plt.lines!(ax1, D_ps, V_0, color = colors[i], label = labels[i], linestyle = :dash)
+        Plt.lines!(
+            ax1,
+            D_ps,
+            V_0,
+            color = colors[i],
+            label = labels[i],
+            linestyle = :dash,
+        )
     end
 
     # F_r = 1
-    ax2 =
-        Plt.Axis(fig[1:7, 10:18], title = "Fᵣ = 1", xlabel = "Dₚ (m)", ylabel = "V (m s⁻¹)")
+    ax2 = Plt.Axis(
+        fig[1:7, 10:18],
+        title = "Fᵣ = 1",
+        xlabel = "Dₚ (m)",
+        ylabel = "V (m s⁻¹)",
+    )
     for i in 1:4
         D_ps, V_1, V_1_ϕ = get_values(
             p3,
@@ -104,7 +121,14 @@ function fig1()
             res,
         )
         Plt.lines!(ax2, D_ps, V_1_ϕ, color = colors[i], label = labels[i])
-        Plt.lines!(ax2, D_ps, V_1, color = colors[i], label = labels[i], linestyle = :dash)
+        Plt.lines!(
+            ax2,
+            D_ps,
+            V_1,
+            color = colors[i],
+            label = labels[i],
+            linestyle = :dash,
+        )
     end
 
     Plt.Legend(
diff --git a/docs/src/plots/CloudDiagnostics.jl b/docs/src/plots/CloudDiagnostics.jl
index 6f3228530..f952335f0 100644
--- a/docs/src/plots/CloudDiagnostics.jl
+++ b/docs/src/plots/CloudDiagnostics.jl
@@ -15,8 +15,13 @@ rain = CMP.Rain(FT)
 cloud_liquid = CMP.CloudLiquid(FT)
 cloud_ice = CMP.CloudIce(FT)
 
-override_file =
-    joinpath(pkgdir(CloudMicrophysics), "src", "parameters", "toml", "SB2006_limiters.toml")
+override_file = joinpath(
+    pkgdir(CloudMicrophysics),
+    "src",
+    "parameters",
+    "toml",
+    "SB2006_limiters.toml",
+)
 toml_dict = CP.create_toml_dict(FT; override_file)
 SB = CMP.SB2006(toml_dict)
 SB_no_limiters = CMP.SB2006(toml_dict, false)
diff --git a/docs/src/plots/Frostenberg_fig1.jl b/docs/src/plots/Frostenberg_fig1.jl
index 9b8c5be28..b1aa14baa 100644
--- a/docs/src/plots/Frostenberg_fig1.jl
+++ b/docs/src/plots/Frostenberg_fig1.jl
@@ -13,8 +13,8 @@ INPC_range = 10.0 .^ (range(-5, stop = 7, length = 500)) #ice nucleating particl
 
 frequency = [
     IN.INP_concentration_frequency(ip, INPC, T) > 0.015 ?
-    IN.INP_concentration_frequency(ip, INPC, T) : missing for INPC in INPC_range,
-    T in T_range
+    IN.INP_concentration_frequency(ip, INPC, T) : missing for
+    INPC in INPC_range, T in T_range
 ]
 mu = [exp(IN.INP_concentration_mean(T)) for T in T_range]
 
@@ -33,7 +33,13 @@ PL.contourf(
     lw = 0,
 )
 
-PL.plot!(T_range, mu, label = "median INPC", legend = :topright, color = :darkred)
+PL.plot!(
+    T_range,
+    mu,
+    label = "median INPC",
+    legend = :topright,
+    color = :darkred,
+)
 
 PL.plot!(
     repeat([257], 50),
diff --git a/docs/src/plots/HomFreezingPlots.jl b/docs/src/plots/HomFreezingPlots.jl
index bdf10ad7f..aaac07e47 100644
--- a/docs/src/plots/HomFreezingPlots.jl
+++ b/docs/src/plots/HomFreezingPlots.jl
@@ -18,9 +18,18 @@ T_range = range(229.0, stop = 234.5, length = 50)  # air temperature
 x_sulph = Vector{FT}([0.03, 0.04, 0.06])           # wt% sulphuric acid in droplets
 
 # Solving for Δa and J values
-Δa1 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[1], T) - CMO.a_w_ice(tps, T) for T in T_range]
-Δa2 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[2], T) - CMO.a_w_ice(tps, T) for T in T_range]
-Δa3 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[3], T) - CMO.a_w_ice(tps, T) for T in T_range]
+Δa1 = [
+    CMO.a_w_xT(H2SO4_prs, tps, x_sulph[1], T) - CMO.a_w_ice(tps, T) for
+    T in T_range
+]
+Δa2 = [
+    CMO.a_w_xT(H2SO4_prs, tps, x_sulph[2], T) - CMO.a_w_ice(tps, T) for
+    T in T_range
+]
+Δa3 = [
+    CMO.a_w_xT(H2SO4_prs, tps, x_sulph[3], T) - CMO.a_w_ice(tps, T) for
+    T in T_range
+]
 
 J1 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa1)
 J2 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa2)
diff --git a/docs/src/plots/KnopfAlpert2013_fig5.jl b/docs/src/plots/KnopfAlpert2013_fig5.jl
index 57ef9419c..00d3550db 100644
--- a/docs/src/plots/KnopfAlpert2013_fig5.jl
+++ b/docs/src/plots/KnopfAlpert2013_fig5.jl
@@ -14,8 +14,16 @@ illite = CMP.Illite(FT)    # dust type
 
 # Knopf and Alpert 2013 Figure 5A: Cirrus
 # https://doi.org/10.1039/C3FD00035D
-temp =
-    [228.20357, 228.33571, 228.50357, 228.75000, 228.92143, 229.16786, 229.39643, 229.52143]
+temp = [
+    228.20357,
+    228.33571,
+    228.50357,
+    228.75000,
+    228.92143,
+    229.16786,
+    229.39643,
+    229.52143,
+]
 KA13_fig5A_J = [
     70170382.86704,
     12101528.74384,
@@ -43,7 +51,12 @@ J_ABIFM = @. CMI.ABIFM_J(illite, Δa_w) * 1e-4 # converted from SI units to cm^-
 
 # Plot results
 fig = MK.Figure(resolution = (800, 600))
-ax1 = MK.Axis(fig[1, 1], ylabel = "J_het [cm^-2 s^-1]", xlabel = "Temp [K]", yscale = log10)
+ax1 = MK.Axis(
+    fig[1, 1],
+    ylabel = "J_het [cm^-2 s^-1]",
+    xlabel = "Temp [K]",
+    yscale = log10,
+)
 
 MK.ylims!(ax1, FT(1), FT(1e8))
 MK.lines!(
diff --git a/docs/src/plots/MarshallPalmer_distribution.jl b/docs/src/plots/MarshallPalmer_distribution.jl
index 27c917192..078d19690 100644
--- a/docs/src/plots/MarshallPalmer_distribution.jl
+++ b/docs/src/plots/MarshallPalmer_distribution.jl
@@ -36,10 +36,14 @@ qᵣ_3 = FT(1.0e-3)
 
 r_range = range(25e-6, stop = 1e-2, length = 1000)
 
-n_r_0 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_0, ρ, r) for r in r_range]
-n_r_1 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_1, ρ, r) for r in r_range]
-n_r_2 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_2, ρ, r) for r in r_range]
-n_r_3 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]
+n_r_0 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_0, ρ, r) for r in r_range]
+n_r_1 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_1, ρ, r) for r in r_range]
+n_r_2 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_2, ρ, r) for r in r_range]
+n_r_3 =
+    [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]
 
 fig = MK.Figure(resolution = (1100, 600))
 ax1 = MK.Axis(
diff --git a/docs/src/plots/Microphysics1M_plots.jl b/docs/src/plots/Microphysics1M_plots.jl
index a4975661d..73fcaf4cf 100644
--- a/docs/src/plots/Microphysics1M_plots.jl
+++ b/docs/src/plots/Microphysics1M_plots.jl
@@ -171,8 +171,15 @@ q_ice = 1e-6
 PL.plot(
     q_rain_range * 1e3,
     [
-        CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for
-        q_rai in q_rain_range
+        CM1.accretion_rain_sink(
+            rain,
+            ice,
+            Blk1MVel.rain,
+            ce,
+            q_ice,
+            q_rai,
+            ρ_air,
+        ) for q_rai in q_rain_range
     ],
     linewidth = 3,
     xlabel = "q_rain or q_snow [g/kg]",
@@ -183,8 +190,15 @@ q_ice = 1e-5
 PL.plot!(
     q_rain_range * 1e3,
     [
-        CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for
-        q_rai in q_rain_range
+        CM1.accretion_rain_sink(
+            rain,
+            ice,
+            Blk1MVel.rain,
+            ce,
+            q_ice,
+            q_rai,
+            ρ_air,
+        ) for q_rai in q_rain_range
     ],
     linewidth = 3,
     xlabel = "q_rain or q_snow [g/kg]",
@@ -195,8 +209,15 @@ q_ice = 1e-4
 PL.plot!(
     q_rain_range * 1e3,
     [
-        CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for
-        q_rai in q_rain_range
+        CM1.accretion_rain_sink(
+            rain,
+            ice,
+            Blk1MVel.rain,
+            ce,
+            q_ice,
+            q_rai,
+            ρ_air,
+        ) for q_rai in q_rain_range
     ],
     linewidth = 3,
     xlabel = "q_rain or q_snow [g/kg]",
@@ -338,8 +359,16 @@ R = TD.gas_constant_air(tps, q)
 PL.plot(
     q_rain_range * 1e3,
     [
-        CM1.evaporation_sublimation(rain, Blk1MVel.rain, aps, tps, q, q_rai, ρ, T) for
-        q_rai in q_rain_range
+        CM1.evaporation_sublimation(
+            rain,
+            Blk1MVel.rain,
+            aps,
+            tps,
+            q,
+            q_rai,
+            ρ,
+            T,
+        ) for q_rai in q_rain_range
     ],
     xlabel = "q_rain [g/kg]",
     linewidth = 3,
@@ -371,8 +400,16 @@ R = TD.gas_constant_air(tps, q)
 PL.plot(
     q_snow_range * 1e3,
     [
-        CM1.evaporation_sublimation(snow, Blk1MVel.snow, aps, tps, q, q_sno, ρ, T) for
-        q_sno in q_snow_range
+        CM1.evaporation_sublimation(
+            snow,
+            Blk1MVel.snow,
+            aps,
+            tps,
+            q,
+            q_sno,
+            ρ,
+            T,
+        ) for q_sno in q_snow_range
     ],
     xlabel = "q_snow [g/kg]",
     linewidth = 3,
@@ -396,8 +433,16 @@ R = TD.gas_constant_air(tps, q)
 PL.plot!(
     q_snow_range * 1e3,
     [
-        CM1.evaporation_sublimation(snow, Blk1MVel.snow, aps, tps, q, q_sno, ρ, T) for
-        q_sno in q_snow_range
+        CM1.evaporation_sublimation(
+            snow,
+            Blk1MVel.snow,
+            aps,
+            tps,
+            q,
+            q_sno,
+            ρ,
+            T,
+        ) for q_sno in q_snow_range
     ],
     xlabel = "q_snow [g/kg]",
     linewidth = 3,
diff --git a/docs/src/plots/Microphysics2M_plots.jl b/docs/src/plots/Microphysics2M_plots.jl
index b8b69633a..1cace8aa7 100644
--- a/docs/src/plots/Microphysics2M_plots.jl
+++ b/docs/src/plots/Microphysics2M_plots.jl
@@ -26,7 +26,9 @@ const liquid = CMP.CloudLiquid(FT)
 const rain = CMP.Rain(FT)
 const blk1mvel = CMP.Blk1MVelType(FT)
 
-include(joinpath(pkgdir(CloudMicrophysics), "docs", "src", "plots", "Wooddata.jl"))
+include(
+    joinpath(pkgdir(CloudMicrophysics), "docs", "src", "plots", "Wooddata.jl"),
+)
 
 # Example values
 q_liq_range = range(1e-8, stop = 1e-3, length = 1000)
@@ -37,44 +39,77 @@ q_rai = 5e-4
 N_d = 1e8
 ρ_air = 1.0 # kg m^-3
 
-q_liq_KK2000 = [CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_B1994 = [CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_TC1980 = [CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_LD2004 = [CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_VarTimeScaleAcnv =
-    [CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
+q_liq_KK2000 = [
+    CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d) for q_liq in q_liq_range
+]
+q_liq_B1994 =
+    [CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d) for q_liq in q_liq_range]
+q_liq_TC1980 = [
+    CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d) for q_liq in q_liq_range
+]
+q_liq_LD2004 = [
+    CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d) for q_liq in q_liq_range
+]
+q_liq_VarTimeScaleAcnv = [
+    CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d) for q_liq in q_liq_range
+]
 q_liq_SB2006 = [
-    CM2.autoconversion(SB2006.acnv, SB2006.pdf_c, q_liq, q_rai, ρ_air, N_d).dq_rai_dt
-    for q_liq in q_liq_range
+    CM2.autoconversion(
+        SB2006.acnv,
+        SB2006.pdf_c,
+        q_liq,
+        q_rai,
+        ρ_air,
+        N_d,
+    ).dq_rai_dt for q_liq in q_liq_range
 ]
-q_liq_K1969 = [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range]
-
-N_d_KK2000 = [CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d) for N_d in N_d_range]
-N_d_B1994 = [CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d) for N_d in N_d_range]
-N_d_TC1980 = [CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d) for N_d in N_d_range]
-N_d_LD2004 = [CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d) for N_d in N_d_range]
+q_liq_K1969 =
+    [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range]
+
+N_d_KK2000 =
+    [CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d) for N_d in N_d_range]
+N_d_B1994 =
+    [CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d) for N_d in N_d_range]
+N_d_TC1980 =
+    [CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d) for N_d in N_d_range]
+N_d_LD2004 =
+    [CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d) for N_d in N_d_range]
 N_d_VarTimeScaleAcnv =
     [CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d) for N_d in N_d_range]
 N_d_SB2006 = [
-    CM2.autoconversion(SB2006.acnv, SB2006.pdf_c, q_liq, q_rai, ρ_air, N_d).dq_rai_dt
-    for N_d in N_d_range
+    CM2.autoconversion(
+        SB2006.acnv,
+        SB2006.pdf_c,
+        q_liq,
+        q_rai,
+        ρ_air,
+        N_d,
+    ).dq_rai_dt for N_d in N_d_range
 ]
 
-accKK2000_q_liq = [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_liq in q_liq_range]
-accB1994_q_liq = [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_liq in q_liq_range]
+accKK2000_q_liq =
+    [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_liq in q_liq_range]
+accB1994_q_liq =
+    [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_liq in q_liq_range]
 accTC1980_q_liq = [CM2.accretion(TC1980, q_liq, q_rai) for q_liq in q_liq_range]
-accSB2006_q_liq =
-    [CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for q_liq in q_liq_range]
+accSB2006_q_liq = [
+    CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for
+    q_liq in q_liq_range
+]
 accK1969_q_liq = [
     CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for
     q_liq in q_liq_range
 ]
 
-accKK2000_q_rai = [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_rai in q_rai_range]
-accB1994_q_rai = [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_rai in q_rai_range]
+accKK2000_q_rai =
+    [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_rai in q_rai_range]
+accB1994_q_rai =
+    [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_rai in q_rai_range]
 accTC1980_q_rai = [CM2.accretion(TC1980, q_liq, q_rai) for q_rai in q_rai_range]
-accSB2006_q_rai =
-    [CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for q_rai in q_rai_range]
+accSB2006_q_rai = [
+    CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for
+    q_rai in q_rai_range
+]
 accK1969_q_rai = [
     CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for
     q_rai in q_rai_range
@@ -97,10 +132,24 @@ l2 = lines!(ax1, q_liq_range * 1e3, q_liq_B1994, color = :green)
 l3 = lines!(ax1, q_liq_range * 1e3, q_liq_TC1980, color = :blue)
 l4 = lines!(ax1, q_liq_range * 1e3, q_liq_LD2004, color = :purple)
 l5 = lines!(ax1, q_liq_range * 1e3, q_liq_K1969, color = :black)
-l6 = lines!(ax1, KK2000_x_q_liq, KK2000_y_q_liq, color = :red, linestyle = :dash)
-l7 = lines!(ax1, B1994_x_q_liq, B1994_y_q_liq, color = :green, linestyle = :dash)
-l8 = lines!(ax1, TC1980_x_q_liq, TC1980_y_q_liq, color = :blue, linestyle = :dash)
-l9 = lines!(ax1, LD2004_x_q_liq, LD2004_y_q_liq, color = :purple, linestyle = :dash)
+l6 =
+    lines!(ax1, KK2000_x_q_liq, KK2000_y_q_liq, color = :red, linestyle = :dash)
+l7 =
+    lines!(ax1, B1994_x_q_liq, B1994_y_q_liq, color = :green, linestyle = :dash)
+l8 = lines!(
+    ax1,
+    TC1980_x_q_liq,
+    TC1980_y_q_liq,
+    color = :blue,
+    linestyle = :dash,
+)
+l9 = lines!(
+    ax1,
+    LD2004_x_q_liq,
+    LD2004_y_q_liq,
+    color = :purple,
+    linestyle = :dash,
+)
 
 l10 = lines!(ax2, N_d_range * 1e-6, N_d_KK2000, color = :red)
 l11 = lines!(ax2, N_d_range * 1e-6, N_d_B1994, color = :green)
@@ -109,7 +158,8 @@ l13 = lines!(ax2, N_d_range * 1e-6, N_d_LD2004, color = :purple)
 l14 = lines!(ax2, KK2000_x_N_d, KK2000_y_N_d, color = :red, linestyle = :dash)
 l15 = lines!(ax2, B1994_x_N_d, B1994_y_N_d, color = :green, linestyle = :dash)
 l16 = lines!(ax2, TC1980_x_N_d, TC1980_y_N_d, color = :blue, linestyle = :dash)
-l17 = lines!(ax2, LD2004_x_N_d, LD2004_y_N_d, color = :purple, linestyle = :dash)
+l17 =
+    lines!(ax2, LD2004_x_N_d, LD2004_y_N_d, color = :purple, linestyle = :dash)
 
 l18 = lines!(ax3, q_liq_range * 1e3, accKK2000_q_liq, color = :red)
 l19 = lines!(ax3, q_liq_range * 1e3, accB1994_q_liq, color = :green)
diff --git a/docs/src/plots/P3LambdaErrorPlots.jl b/docs/src/plots/P3LambdaErrorPlots.jl
index a100a8573..a692b3821 100644
--- a/docs/src/plots/P3LambdaErrorPlots.jl
+++ b/docs/src/plots/P3LambdaErrorPlots.jl
@@ -24,7 +24,8 @@ function λ_diff(F_rim::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3) where {FT}
     # Compute mass density based on input shape parameters
     L_calc = N * P3.L_over_N_gamma(p3, F_rim, F_liq, x, μ, th)
 
-    (λ_calculated,) = P3.distribution_parameter_solver(p3, L_calc, N, ρ_r, F_rim, F_liq)
+    (λ_calculated,) =
+        P3.distribution_parameter_solver(p3, L_calc, N, ρ_r, F_rim, F_liq)
     return abs(λ_ex - λ_calculated)
 end
 
@@ -154,7 +155,14 @@ function μ_approximation_effects(F_rim::FT, ρ_r::FT) where {FT}
         L_over_N = P3.L_over_N_gamma(p3, F_rim, F_liq, log(λs[i]), μs[i], th)
         Ls[i] = L_over_N
         N = FT(1e6)
-        (L, N) = P3.distribution_parameter_solver(p3, L_over_N * N, N, ρ_r, F_rim, F_liq)
+        (L, N) = P3.distribution_parameter_solver(
+            p3,
+            L_over_N * N,
+            N,
+            ρ_r,
+            F_rim,
+            F_liq,
+        )
         λ_solved[i] = L
         μs_approx[i] = P3.DSD_μ_approx(p3, N * L_over_N, N, ρ_r, F_rim, F_liq)
     end
diff --git a/docs/src/plots/P3ShapeSolverPlots.jl b/docs/src/plots/P3ShapeSolverPlots.jl
index 55c2ef81f..ab81c87ab 100644
--- a/docs/src/plots/P3ShapeSolverPlots.jl
+++ b/docs/src/plots/P3ShapeSolverPlots.jl
@@ -32,7 +32,12 @@ function lambda_guess_plot()
                 f[i, j],
                 xlabel = "log(L/N)",
                 ylabel = "log(λ)",
-                title = string("λ vs L/N for F_rim = ", F_rim, " and ρ_r = ", ρ_r),
+                title = string(
+                    "λ vs L/N for F_rim = ",
+                    F_rim,
+                    " and ρ_r = ",
+                    ρ_r,
+                ),
                 height = 300,
                 width = 400,
             )
@@ -41,7 +46,14 @@ function lambda_guess_plot()
             λs = [10^logλ for logλ in logλs]
             th = P3.thresholds(p3, ρ_r, F_rim)
             Ls_over_N = [
-                P3.L_over_N_gamma(p3, F_rim, F_liq, log(λ), P3.DSD_μ(p3, λ), th) for λ in λs
+                P3.L_over_N_gamma(
+                    p3,
+                    F_rim,
+                    F_liq,
+                    log(λ),
+                    P3.DSD_μ(p3, λ),
+                    th,
+                ) for λ in λs
             ]
             guesses = [FT(0) for λ in λs]
 
diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl
index 7c94c1951..cad714dd4 100644
--- a/docs/src/plots/P3TerminalVelocityPlots.jl
+++ b/docs/src/plots/P3TerminalVelocityPlots.jl
@@ -136,8 +136,13 @@ function figure_2()
 
     # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015
 
-    contour_args =
-        (color = :black, labels = true, levels = 3, linewidth = 1.5, labelsize = 18)
+    contour_args = (
+        color = :black,
+        labels = true,
+        levels = 3,
+        linewidth = 1.5,
+        labelsize = 18,
+    )
 
     ax1 = make_axis(fig, 1, 1, "Particle regimes with small Dₘ")
     hm = Plt.contourf!(
diff --git a/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
index 358fad4e1..55f2d62d8 100644
--- a/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
+++ b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
@@ -117,7 +117,13 @@ function fig1()
     (F_riml, F_liql, V_ml, D_ml, D_ml_regimes) =
         get_values(p3, Chen2022, L_l, N_l, ρ_a, xres, yres)
 
-    args = (color = :black, labels = true, levels = 3, linewidth = 1.5, labelsize = 18)
+    args = (
+        color = :black,
+        labels = true,
+        levels = 3,
+        linewidth = 1.5,
+        labelsize = 18,
+    )
 
 
     fig = Plt.Figure()
diff --git a/docs/src/plots/RainEvapoartionSB2006.jl b/docs/src/plots/RainEvapoartionSB2006.jl
index d1ef21711..a897a7b60 100644
--- a/docs/src/plots/RainEvapoartionSB2006.jl
+++ b/docs/src/plots/RainEvapoartionSB2006.jl
@@ -36,10 +36,12 @@ function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T)
         t_star = (FT(6) * x_star / xr)^FT(1 / 3)
         a_vent_0 = av * FT(SF.gamma(-1, t_star)) / FT(6)^FT(-2 / 3)
         b_vent_0 =
-            bv * FT(SF.gamma((-1 / 2) + (3 / 2) * β, t_star)) / FT(6)^FT(β / 2 - 1 / 2)
+            bv * FT(SF.gamma((-1 / 2) + (3 / 2) * β, t_star)) /
+            FT(6)^FT(β / 2 - 1 / 2)
 
         a_vent_1 = av * SF.gamma(FT(2)) / FT(6)^FT(1 / 3)
-        b_vent_1 = bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
+        b_vent_1 =
+            bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
 
         N_Re = α * xr^β * sqrt(ρ0 / ρ) * Dr / ν_air
         Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
diff --git a/docs/src/plots/Riccobono_mixed_nucleation_plots.jl b/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
index b2dc04afb..e97df0b6a 100644
--- a/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
+++ b/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
@@ -19,7 +19,11 @@ nucleation_rates = map(bioOxOrg_concentrations) do bioOxOrg_conc
 end
 
 Plots.plot(xaxis = :log, yaxis = :log, lw = 3)
-Plots.plot!(bioOxOrg_concentrations, nucleation_rates, label = "Riccobono parameterization")
+Plots.plot!(
+    bioOxOrg_concentrations,
+    nucleation_rates,
+    label = "Riccobono parameterization",
+)
 CLOUD_points = [
     (5360344.038850841, 0.09948622490099472),
     (8968652.585232794, 0.3074108277542233),
diff --git a/docs/src/plots/TerminalVelocity2M.jl b/docs/src/plots/TerminalVelocity2M.jl
index 34a8401fc..c2cc6b602 100644
--- a/docs/src/plots/TerminalVelocity2M.jl
+++ b/docs/src/plots/TerminalVelocity2M.jl
@@ -17,7 +17,8 @@ const SB2006Vel = CMP.SB2006VelType(FT)
 q_rain_random = rand(10000) .* 5e-3
 N_rain_random = 10.0 .^ (1.0 .+ 6.0 .* rand(10000))
 r_mean =
-    ((ρ_air .* q_rain_random ./ N_rain_random) / 1000.0 * 3 / 4 / pi) .^ (1.0 / 3) .* 1e6
+    ((ρ_air .* q_rain_random ./ N_rain_random) / 1000.0 * 3 / 4 / pi) .^
+    (1.0 / 3) .* 1e6
 
 #! format: off
 
diff --git a/docs/src/plots/Thersholds_transitions.jl b/docs/src/plots/Thersholds_transitions.jl
index c758dd96d..0992a78fe 100644
--- a/docs/src/plots/Thersholds_transitions.jl
+++ b/docs/src/plots/Thersholds_transitions.jl
@@ -39,22 +39,36 @@ N_d_range = range(1e7, stop = 1e9, length = 1000)
 ρ_air = 1.0 # kg m^-3
 N_d = 1e8
 
-q_liq_K1969 = [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq) for q_liq in q_liq_range]
-q_liq_K1969_s =
-    [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq, true) for q_liq in q_liq_range]
+q_liq_K1969 =
+    [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq) for q_liq in q_liq_range]
+q_liq_K1969_s = [
+    CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq, true) for
+    q_liq in q_liq_range
+]
 
-q_liq_TC1980 =
-    [CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_TC1980_s =
-    [CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, N_d, true) for q_liq in q_liq_range]
-q_liq_LD2004 =
-    [CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, N_d) for q_liq in q_liq_range]
-q_liq_LD2004_s =
-    [CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, N_d, true) for q_liq in q_liq_range]
+q_liq_TC1980 = [
+    CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, N_d) for
+    q_liq in q_liq_range
+]
+q_liq_TC1980_s = [
+    CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, N_d, true) for
+    q_liq in q_liq_range
+]
+q_liq_LD2004 = [
+    CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, N_d) for
+    q_liq in q_liq_range
+]
+q_liq_LD2004_s = [
+    CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, N_d, true) for
+    q_liq in q_liq_range
+]
 
-N_d_B1994 = [CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d) for N_d in N_d_range]
-N_d_B1994_s =
-    [CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d, true) for N_d in N_d_range]
+N_d_B1994 =
+    [CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d) for N_d in N_d_range]
+N_d_B1994_s = [
+    CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d, true) for
+    N_d in N_d_range
+]
 
 PL.plot(
     q_liq_range * 1e3,
@@ -64,7 +78,12 @@ PL.plot(
     ylabel = "autoconversion rate [1/s]",
     label = "K1969 without smoothing",
 )
-PL.plot!(q_liq_range * 1e3, q_liq_K1969_s, linewidth = 2, label = "K1969 with smoothing")
+PL.plot!(
+    q_liq_range * 1e3,
+    q_liq_K1969_s,
+    linewidth = 2,
+    label = "K1969 with smoothing",
+)
 PL.savefig("q_liq_K1969.svg") # hide
 
 PL.plot(
@@ -77,9 +96,24 @@ PL.plot(
     yaxis = :log,
     ylim = (1e-10, 1e-5),
 )
-PL.plot!(q_liq_range * 1e3, q_liq_TC1980_s, linewidth = 2, label = "TC1980 with smoothing")
-PL.plot!(q_liq_range * 1e3, q_liq_LD2004, linewidth = 2, label = "LD2004 without smoothing")
-PL.plot!(q_liq_range * 1e3, q_liq_LD2004_s, linewidth = 2, label = "LD2004 with smoothing")
+PL.plot!(
+    q_liq_range * 1e3,
+    q_liq_TC1980_s,
+    linewidth = 2,
+    label = "TC1980 with smoothing",
+)
+PL.plot!(
+    q_liq_range * 1e3,
+    q_liq_LD2004,
+    linewidth = 2,
+    label = "LD2004 without smoothing",
+)
+PL.plot!(
+    q_liq_range * 1e3,
+    q_liq_LD2004_s,
+    linewidth = 2,
+    label = "LD2004 with smoothing",
+)
 PL.savefig("q_liq_TC1980_LD2004.svg") # hide
 
 PL.plot(
@@ -93,5 +127,10 @@ PL.plot(
     yaxis = :log,
     ylim = (1e-13, 1e-5),
 )
-PL.plot!(N_d_range * 1e-6, N_d_B1994_s, linewidth = 2, label = "B1994 with smoothing")
+PL.plot!(
+    N_d_range * 1e-6,
+    N_d_B1994_s,
+    linewidth = 2,
+    label = "B1994 with smoothing",
+)
 PL.savefig("N_d_B1994.svg") # hide
diff --git a/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl b/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
index b93af8709..eb14f5dde 100644
--- a/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
+++ b/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
@@ -30,9 +30,9 @@ function vehkamaki_nucleation_timestep(rh, temp, so4)
     x_crit = (
         0.740997 - 0.00266379 * temp - 0.00349998 * ln_so4 +
         0.0000504022 * temp * ln_so4 +
-        0.00201048 * log_rh - 0.000183289 * temp * log_rh + 0.00157407 * log_rh^2 -
-        0.0000179059 * temp * log_rh^2 + 0.000184403 * log_rh^3 -
-        1.50345e-6 * temp * log_rh^3
+        0.00201048 * log_rh - 0.000183289 * temp * log_rh +
+        0.00157407 * log_rh^2 - 0.0000179059 * temp * log_rh^2 +
+        0.000184403 * log_rh^3 - 1.50345e-6 * temp * log_rh^3
     )
 
     a_T = (
@@ -46,24 +46,24 @@ function vehkamaki_nucleation_timestep(rh, temp, so4)
         15.7963 / x_crit
     )
     c_T = (
-        -0.215554 - 0.0810269 * temp + 0.00143581 * temp^2 - 4.7758e-6 * temp^3 -
-        2.91297 / x_crit
+        -0.215554 - 0.0810269 * temp + 0.00143581 * temp^2 -
+        4.7758e-6 * temp^3 - 2.91297 / x_crit
     )
     d_T = (
-        -3.58856 + 0.049508 * temp - 0.00021382 * temp^2 + 3.10801e-7 * temp^3 -
-        0.0293333 / x_crit
+        -3.58856 + 0.049508 * temp - 0.00021382 * temp^2 + 3.10801e-7 * temp^3 - 0.0293333 / x_crit
     )
     e_T = (
-        1.14598 - 0.600796 * temp + 0.00864245 * temp^2 - 0.0000228947 * temp^3 -
-        8.44985 / x_crit
+        1.14598 - 0.600796 * temp + 0.00864245 * temp^2 -
+        0.0000228947 * temp^3 - 8.44985 / x_crit
     )
     f_T = (
-        2.15855 + 0.0808121 * temp - 0.000407382 * temp^2 - 4.01957e-7 * temp^3 +
-        0.721326 / x_crit
+        2.15855 + 0.0808121 * temp - 0.000407382 * temp^2 -
+        4.01957e-7 * temp^3 + 0.721326 / x_crit
     )
     g_T = (
-        1.6241 - 0.0160106 * temp + 0.0000377124 * temp^2 + 3.21794e-8 * temp^3 -
-        0.0113255 / x_crit
+        1.6241 - 0.0160106 * temp +
+        0.0000377124 * temp^2 +
+        3.21794e-8 * temp^3 - 0.0113255 / x_crit
     )
     h_T = (
         9.71682 - 0.115048 * temp +
@@ -72,12 +72,12 @@ function vehkamaki_nucleation_timestep(rh, temp, so4)
         0.71186 / x_crit
     )
     i_T = (
-        -1.05611 + 0.00903378 * temp - 0.0000198417 * temp^2 + 2.46048e-8 * temp^3 -
-        0.0579087 / x_crit
+        -1.05611 + 0.00903378 * temp - 0.0000198417 * temp^2 +
+        2.46048e-8 * temp^3 - 0.0579087 / x_crit
     )
     j_T = (
-        -0.148712 + 0.00283508 * temp - 9.24619e-6 * temp^2 + 5.00427e-9 * temp^3 -
-        0.0127081 / x_crit
+        -0.148712 + 0.00283508 * temp - 9.24619e-6 * temp^2 +
+        5.00427e-9 * temp^3 - 0.0127081 / x_crit
     )
 
     # Nucleation rate 1/(m^3 s)
diff --git a/docs/src/plots/linear_HOM_J.jl b/docs/src/plots/linear_HOM_J.jl
index a7ecc602c..997c272c2 100644
--- a/docs/src/plots/linear_HOM_J.jl
+++ b/docs/src/plots/linear_HOM_J.jl
@@ -14,7 +14,8 @@ log10KoopJ = @. log10(KoopJ .* 1e-6)
 # Linear fit on Koop 2000
 M = [ones(length(Koop_Δa)) Koop_Δa]
 linear_coeffs = M \ log10KoopJ
-new_log10J = [linear_coeffs[2] * Delta_a + linear_coeffs[1] for Delta_a in Koop_Δa]
+new_log10J =
+    [linear_coeffs[2] * Delta_a + linear_coeffs[1] for Delta_a in Koop_Δa]
 
 # Plotting J vs Δa
 fig = MK.Figure(resolution = (800, 600), fontsize = 18)
diff --git a/ext/Common.jl b/ext/Common.jl
index c5642950e..893f00139 100644
--- a/ext/Common.jl
+++ b/ext/Common.jl
@@ -27,7 +27,10 @@ function get_num_modes(data_row::NamedTuple)
     end
 end
 
-function read_aerosol_dataset(dataset_filename::String, Y_name = :mode_1_act_frac_S_interp)
+function read_aerosol_dataset(
+    dataset_filename::String,
+    Y_name = :mode_1_act_frac_S_interp,
+)
     initial_data = DF.DataFrame(CSV.File(dataset_filename))
     df = filter(row -> row.S_max > 0 && row.S_max < 0.2, initial_data)
     selected_columns_X = []
@@ -43,7 +46,10 @@ function read_aerosol_dataset(dataset_filename::String, Y_name = :mode_1_act_fra
             ]),
         )
     end
-    append!(selected_columns_X, [:velocity, :initial_temperature, :initial_pressure])
+    append!(
+        selected_columns_X,
+        [:velocity, :initial_temperature, :initial_pressure],
+    )
     X = df[:, selected_columns_X]
     Y = df[:, Y_name]
     return (X, Y, initial_data)
@@ -52,10 +58,14 @@ end
 function preprocess_aerosol_data(X::DataFrame)
     num_modes = get_num_modes(X)
     for i in 1:num_modes
-        X = DF.transform(X, Symbol("mode_$(i)_N") => ByRow(log) => Symbol("mode_$(i)_N"))
         X = DF.transform(
             X,
-            Symbol("mode_$(i)_mean") => ByRow(log) => Symbol("mode_$(i)_mean"),
+            Symbol("mode_$(i)_N") => ByRow(log) => Symbol("mode_$(i)_N"),
+        )
+        X = DF.transform(
+            X,
+            Symbol("mode_$(i)_mean") =>
+                ByRow(log) => Symbol("mode_$(i)_mean"),
         )
     end
     X = DF.transform(X, :velocity => ByRow(log) => :velocity)
@@ -103,8 +113,9 @@ function get_ARG_act_frac(
     q_vap = vapor_mix_ratio / (vapor_mix_ratio + 1)
     q = TD.PhasePartition(FT(q_vap), FT(0), FT(0))
 
-    return collect(AA.N_activated_per_mode(ap, ad, aip, tps, FT(T), FT(p), FT(w), q)) ./
-           mode_Ns
+    return collect(
+        AA.N_activated_per_mode(ap, ad, aip, tps, FT(T), FT(p), FT(w), q),
+    ) ./ mode_Ns
 end
 
 function get_ARG_act_frac(
@@ -114,7 +125,9 @@ function get_ARG_act_frac(
     tps::TDP.ThermodynamicsParameters,
     FT::DataType,
 )
-    return transpose(hcat(get_ARG_act_frac.(NamedTuple.(eachrow(X)), ap, aip, tps, FT)...))
+    return transpose(
+        hcat(get_ARG_act_frac.(NamedTuple.(eachrow(X)), ap, aip, tps, FT)...),
+    )
 end
 
 function preprocess_aerosol_data_with_ARG_act_frac(
@@ -129,7 +142,8 @@ function preprocess_aerosol_data_with_ARG_act_frac(
     X = DF.transform(
         X,
         AsTable(All()) =>
-            ByRow(x -> f(x)) => [Symbol("mode_$(i)_ARG_act_frac") for i in 1:num_modes],
+            ByRow(x -> f(x)) =>
+                [Symbol("mode_$(i)_ARG_act_frac") for i in 1:num_modes],
     )
     return preprocess_aerosol_data(X)
 end
diff --git a/ext/EmulatorModelsExt.jl b/ext/EmulatorModelsExt.jl
index 21719afb5..e5824c502 100644
--- a/ext/EmulatorModelsExt.jl
+++ b/ext/EmulatorModelsExt.jl
@@ -51,8 +51,11 @@ function AA.N_activated_per_mode(
                 Symbol("mode_$(j)_kappa") => hygro[modes_perm[j]],
             ) for j in 1:AM.n_modes(ad)
         ]
-        additional_data =
-            (; :velocity => w, :initial_temperature => T, :initial_pressure => p)
+        additional_data = (;
+            :velocity => w,
+            :initial_temperature => T,
+            :initial_pressure => p,
+        )
         X = DF.DataFrame([merge(reduce(merge, per_mode_data), additional_data)])
         max(FT(0), min(FT(1), MLJ.predict(machine, X)[1])) * ad.modes[i].N
     end
@@ -93,7 +96,9 @@ end
     - `ekp_params` - parameters from the trained Ensemble Kalman Process
 Returns a calibrated set of aerosol activation parameters
 """
-function CMP.AerosolActivationParameters(ekp_params::Array{FT}) where {FT <: Real}
+function CMP.AerosolActivationParameters(
+    ekp_params::Array{FT},
+) where {FT <: Real}
     default_param_set = CMP.AerosolActivationParameters(FT)
     (f1, f2, g1, g2, p1, p2) = FT.(ekp_params)
     cur_values = (;
diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl
index 80aa36e82..184990973 100644
--- a/papers/ice_nucleation_2024/AIDA_calibrations.jl
+++ b/papers/ice_nucleation_2024/AIDA_calibrations.jl
@@ -38,8 +38,11 @@ function moving_average(data, n)
 end
 
 # Defining data names, start/end times, etc.
-data_file_names =
-    [["in05_17_aida.edf", "in05_18_aida.edf"], ["in07_01_aida.edf"], ["in07_19_aida.edf"]]
+data_file_names = [
+    ["in05_17_aida.edf", "in05_18_aida.edf"],
+    ["in07_01_aida.edf"],
+    ["in07_19_aida.edf"],
+]
 batch_names = ["IN05", "IN0701", "IN0719"]
 end_sim = 25               # Loss func looks at last end_sim timesteps only
 start_time_list =          # freezing onset
diff --git a/parcel/Example_Frostenberg_Immersion_Freezing.jl b/parcel/Example_Frostenberg_Immersion_Freezing.jl
index d7d64621d..77035fe88 100644
--- a/parcel/Example_Frostenberg_Immersion_Freezing.jl
+++ b/parcel/Example_Frostenberg_Immersion_Freezing.jl
@@ -146,7 +146,12 @@ ax4 = MK.Axis(fig[2, 2]; ylabel = "q_liq [g/kg]")
 ax5 = MK.Axis(fig[3, 1]; xlabel = "Time [min]", ylabel = "N_ice [1/cm3]")
 ax6 = MK.Axis(fig[3, 2]; xlabel = "Time [min]", ylabel = "N_liq [1/cm3]")
 MK.xlims!.(fig.content, 0, t_max / 60)
-ax7 = MK.Axis(fig[4, 1:2]; xlabel = "T [K]", ylabel = "INPC [1/m3]", yscale = log10)
+ax7 = MK.Axis(
+    fig[4, 1:2];
+    xlabel = "T [K]",
+    ylabel = "INPC [1/m3]",
+    yscale = log10,
+)
 # top axis with time
 ax7_top = MK.Axis(fig[4, 1:2]; xaxisposition = :top, xlabel = "Time [min]")
 MK.hidespines!(ax7_top)
@@ -157,12 +162,25 @@ colors = [:red, :green, :orange, :limegreen]
 colors_mean = [:black, :gray]
 
 function plot_results(sol, t, ll, cl, ls)
-    MK.lines!(ax1, t / 60, S_i.(tps, sol[3, :], sol[1, :]) .- 1, linestyle = ls, color = cl)
+    MK.lines!(
+        ax1,
+        t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        linestyle = ls,
+        color = cl,
+    )
     MK.lines!(ax2, t / 60, sol[3, :], linestyle = ls, color = cl)
     MK.lines!(ax3, t / 60, sol[6, :] * 1e3, linestyle = ls, color = cl)
     MK.lines!(ax4, t / 60, sol[5, :] * 1e3, linestyle = ls, color = cl)
     MK.lines!(ax5, t / 60, sol[9, :] * 1e-6, linestyle = ls, color = cl)
-    MK.lines!(ax6, t / 60, sol[8, :] * 1e-6, linestyle = ls, color = cl, label = ll)
+    MK.lines!(
+        ax6,
+        t / 60,
+        sol[8, :] * 1e-6,
+        linestyle = ls,
+        color = cl,
+        label = ll,
+    )
 end
 for (sol, t, ll, cl) in zip(results_stoch, time_stoch, labels_stoch, colors)
     ls = :solid
diff --git a/parcel/Example_Immersion_Freezing.jl b/parcel/Example_Immersion_Freezing.jl
index 41e019a5b..05afa232b 100644
--- a/parcel/Example_Immersion_Freezing.jl
+++ b/parcel/Example_Immersion_Freezing.jl
@@ -71,7 +71,12 @@ for DSD in size_distribution_list
     local sol = run_parcel(IC, FT(0), t_max, params)
 
     # Plot results
-    MK.lines!(ax1, sol.t / 60, S_i.(tps, sol[3, :], sol[1, :]) .- 1, label = DSD)
+    MK.lines!(
+        ax1,
+        sol.t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        label = DSD,
+    )
     MK.lines!(ax2, sol.t / 60, sol[3, :])
     MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
     MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)
diff --git a/parcel/Example_Jensen_et_al_2022.jl b/parcel/Example_Jensen_et_al_2022.jl
index 871be137a..b455c4eee 100644
--- a/parcel/Example_Jensen_et_al_2022.jl
+++ b/parcel/Example_Jensen_et_al_2022.jl
@@ -101,7 +101,14 @@ MK.lines!(
     color = :blue,
     linewidth = 2,
 )
-MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red, linewidth = 2) # liq saturation
+MK.lines!(
+    ax1,
+    sol.t,
+    (sol[1, :]),
+    label = "liquid",
+    color = :red,
+    linewidth = 2,
+) # liq saturation
 MK.lines!(ax2, sol.t, sol[9, :] * 1e-6, label = "CM.jl", linewidth = 2) # ICNC
 
 MK.axislegend(
diff --git a/parcel/Example_NonEq.jl b/parcel/Example_NonEq.jl
index 1731f6655..a8cf2248d 100644
--- a/parcel/Example_NonEq.jl
+++ b/parcel/Example_NonEq.jl
@@ -43,14 +43,18 @@ qᵥ = mv_v / (md_v + mv_v + ml_v + mi_v)
 qₗ = ml_v / (md_v + mv_v + ml_v + mi_v)
 qᵢ = mi_v / (md_v + mv_v + ml_v + mi_v)
 
-override_file =
-    Dict("condensation_evaporation_timescale" => Dict("value" => τ, "type" => "float"))
+override_file = Dict(
+    "condensation_evaporation_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
 
 liquid_toml_dict = CP.create_toml_dict(FT; override_file)
 liquid = CMP.CloudLiquid(liquid_toml_dict)
 
-override_file =
-    Dict("sublimation_deposition_timescale" => Dict("value" => τ, "type" => "float"))
+override_file = Dict(
+    "sublimation_deposition_timescale" =>
+        Dict("value" => τ, "type" => "float"),
+)
 
 ice_toml_dict = CP.create_toml_dict(FT; override_file)
 
diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl
index f9891c405..b8bbf9757 100644
--- a/parcel/ParcelModel.jl
+++ b/parcel/ParcelModel.jl
@@ -265,7 +265,8 @@ function run_parcel(IC, t_0, t_end, pp)
     elseif pp.heterogeneous == "P3_het"
         imm_params = P3_het{FT}(pp.ips, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_random"
-        imm_params = Frostenberg_random{FT}(pp.ip, pp.sampling_interval, pp.const_dt)
+        imm_params =
+            Frostenberg_random{FT}(pp.ip, pp.sampling_interval, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_mean"
         imm_params = Frostenberg_mean{FT}(pp.ip, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_stochastic"
diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl
index 5dc7e6732..37e3984f0 100644
--- a/parcel/ParcelTendencies.jl
+++ b/parcel/ParcelTendencies.jl
@@ -22,7 +22,15 @@ function aerosol_activation(params::AeroAct, state)
     (; T, Sₗ, Nₐ) = state
     FT = eltype(state)
 
-    ad = AM.Mode_κ(r_nuc, aero_σ_g, Nₐ, (FT(1.0),), (FT(1.0),), (aerosol.M,), (aerosol.κ,))
+    ad = AM.Mode_κ(
+        r_nuc,
+        aero_σ_g,
+        Nₐ,
+        (FT(1.0),),
+        (FT(1.0),),
+        (aerosol.M,),
+        (aerosol.κ,),
+    )
     all_ad = AM.AerosolDistribution((ad,))
 
     smax = (Sₗ - 1) < 0 ? FT(0) : (Sₗ - 1)
@@ -48,7 +56,12 @@ function deposition_nucleation(params::MohlerAF, state, dY)
 
     AF =
         Sᵢ >= ips.deposition.Sᵢ_max ? FT(0) :
-        AF = CMI_het.dust_activated_number_fraction(aerosol, ips.deposition, Sᵢ, T)
+        AF = CMI_het.dust_activated_number_fraction(
+            aerosol,
+            ips.deposition,
+            Sᵢ,
+            T,
+        )
     return max(FT(0), AF * Nₐ - Nᵢ) / const_dt
 end
 
@@ -63,7 +76,14 @@ function deposition_nucleation(params::MohlerRate, state, dY)
         @warn("Supersaturation exceeds Sᵢ_max. No dust will be activated.")
         dNi_dt = FT(0)
     else
-        dNi_dt = CMI_het.MohlerDepositionRate(aerosol, ips.deposition, Sᵢ, T, dSᵢdt, Nₐ)
+        dNi_dt = CMI_het.MohlerDepositionRate(
+            aerosol,
+            ips.deposition,
+            Sᵢ,
+            T,
+            dSᵢdt,
+            Nₐ,
+        )
     end
 
     return min(max(dNi_dt, FT(0)), Nₐ / const_dt)
diff --git a/src/AerosolActivation.jl b/src/AerosolActivation.jl
index e75084652..afa0ca9aa 100644
--- a/src/AerosolActivation.jl
+++ b/src/AerosolActivation.jl
@@ -33,7 +33,10 @@ export mean_hygroscopicity_parameter,
 
 Returns a curvature coefficient.
 """
-function coeff_of_curvature(ap::CMP.AerosolActivationParameters, T::FT) where {FT}
+function coeff_of_curvature(
+    ap::CMP.AerosolActivationParameters,
+    T::FT,
+) where {FT}
     return FT(2) * ap.σ * ap.M_w / ap.ρ_w / ap.R / T
 end
 
@@ -166,7 +169,9 @@ function max_supersaturation(
         g::FT = ap.g1 + ap.g2 * log(mode_i.stdev)
         η::FT = (α * w / G)^FT(3 / 2) / (FT(2 * pi) * ap.ρ_w * γ * mode_i.N)
 
-        tmp += 1 / (Sm[i])^2 * (f * (ζ / η)^ap.p1 + g * (Sm[i]^2 / (η + 3 * ζ))^ap.p2)
+        tmp +=
+            1 / (Sm[i])^2 *
+            (f * (ζ / η)^ap.p1 + g * (Sm[i]^2 / (η + 3 * ζ))^ap.p2)
     end
 
     return FT(1) / sqrt(tmp)
@@ -248,7 +253,8 @@ function M_activated_per_mode(
 
         u_i = 2 * log(sm[i] / smax) / 3 / sqrt(2) / log(mode_i.stdev)
 
-        avg_molar_mass_i * 1 / 2 * (1 - SF.erf(u_i - 3 * sqrt(2) / 2 * log(mode_i.stdev)))
+        avg_molar_mass_i * 1 / 2 *
+        (1 - SF.erf(u_i - 3 * sqrt(2) / 2 * log(mode_i.stdev)))
     end
 end
 
diff --git a/src/CloudDiagnostics.jl b/src/CloudDiagnostics.jl
index b118dbbc9..2cd51ee28 100644
--- a/src/CloudDiagnostics.jl
+++ b/src/CloudDiagnostics.jl
@@ -24,7 +24,11 @@ based on the assumed rain particle size distribution.
 Normalized by the reflectivty of 1 millimiter drop in a volume of 1m3.
 The values are clipped at -150 dBZ.
 """
-function radar_reflectivity_1M((; pdf, mass)::CMP.Rain{FT}, q::FT, ρ::FT) where {FT}
+function radar_reflectivity_1M(
+    (; pdf, mass)::CMP.Rain{FT},
+    q::FT,
+    ρ::FT,
+) where {FT}
 
     # change units for accuracy
     n0 = CM1.get_n0(pdf) * FT(1e-12)
@@ -75,11 +79,16 @@ function radar_reflectivity_2M(
 
     Zc =
         Bc < eps(FT) ? FT(0) :
-        Ac * C^(νc + 1) * (Bc * C^μc)^(-(3 + νc) / μc) * SF.gamma((3 + νc) / μc) / μc *
-        FT(10)^(-2 * χ)
+        Ac *
+        C^(νc + 1) *
+        (Bc * C^μc)^(-(3 + νc) / μc) *
+        SF.gamma((3 + νc) / μc) / μc * FT(10)^(-2 * χ)
     Zr =
         Br < eps(FT) ? FT(0) :
-        Ar * C^(νr + 1) * (Br * C^μr)^(-(3 + νr) / μr) * SF.gamma((3 + νr) / μr) / μr
+        Ar *
+        C^(νr + 1) *
+        (Br * C^μr)^(-(3 + νr) / μr) *
+        SF.gamma((3 + νr) / μr) / μr
 
     return max(FT(-150), 10 * (log10(max(FT(0), Zc + Zr)) - log_10_Z₀))
 end
@@ -118,11 +127,16 @@ function effective_radius_2M(
 
     M3_c =
         Bc == 0 ? FT(0) :
-        Ac * C^(νc + 1) * (Bc * C^μc)^(-(2 + νc) / μc) * SF.gamma((2 + νc) / μc) / μc /
-        FT(10)^χ
+        Ac *
+        C^(νc + 1) *
+        (Bc * C^μc)^(-(2 + νc) / μc) *
+        SF.gamma((2 + νc) / μc) / μc / FT(10)^χ
     M3_r =
         Br == 0 ? FT(0) :
-        Ar * C^(νr + 1) * (Br * C^μr)^(-(2 + νr) / μr) * SF.gamma((2 + νr) / μr) / μr
+        Ar *
+        C^(νr + 1) *
+        (Br * C^μr)^(-(2 + νr) / μr) *
+        SF.gamma((2 + νr) / μr) / μr
 
     M2_c =
         Bc == 0 ? FT(0) :
@@ -166,7 +180,10 @@ function effective_radius_Liu_Hallet_97(
     k = FT(0.8)
     r_vol =
         ((N_liq + N_rai) < eps(FT)) ? FT(0) :
-        ((FT(3) * (q_liq + q_rai) * ρ_air) / (FT(4) * π * ρw * (N_liq + N_rai)))^FT(1 / 3)
+        (
+            (FT(3) * (q_liq + q_rai) * ρ_air) /
+            (FT(4) * π * ρw * (N_liq + N_rai))
+        )^FT(1 / 3)
 
     return r_vol / k^FT(1 / 3)
 end
diff --git a/src/Common.jl b/src/Common.jl
index b7cbad5d5..31422821c 100644
--- a/src/Common.jl
+++ b/src/Common.jl
@@ -43,7 +43,8 @@ function G_func(
     L = TD.latent_heat_vapor(tps, T)
     p_vs = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
 
-    return FT(1) / (L / K_therm / T * (L / R_v / T - FT(1)) + R_v * T / D_vapor / p_vs)
+    return FT(1) /
+           (L / K_therm / T * (L / R_v / T - FT(1)) + R_v * T / D_vapor / p_vs)
 end
 function G_func(
     (; K_therm, D_vapor)::CMP.AirProperties{FT},
@@ -55,7 +56,8 @@ function G_func(
     L = TD.latent_heat_sublim(tps, T)
     p_vs = TD.saturation_vapor_pressure(tps, T, TD.Ice())
 
-    return FT(1) / (L / K_therm / T * (L / R_v / T - FT(1)) + R_v * T / D_vapor / p_vs)
+    return FT(1) /
+           (L / K_therm / T * (L / R_v / T - FT(1)) + R_v * T / D_vapor / p_vs)
 end
 
 """
@@ -129,7 +131,18 @@ Returns the saturation vapor pressure above a sulphuric acid solution droplet in
 `x` is, for example, 0.1 if droplets are 10 percent sulphuric acid by weight
 """
 function H2SO4_soln_saturation_vapor_pressure(
-    (; T_max, T_min, w_2, c1, c2, c3, c4, c5, c6, c7)::CMP.H2SO4SolutionParameters{FT},
+    (;
+        T_max,
+        T_min,
+        w_2,
+        c1,
+        c2,
+        c3,
+        c4,
+        c5,
+        c6,
+        c7,
+    )::CMP.H2SO4SolutionParameters{FT},
     x::FT,
     T::FT,
 ) where {FT}
@@ -140,7 +153,8 @@ function H2SO4_soln_saturation_vapor_pressure(
     w_h = w_2 * x
     p_sol =
         exp(
-            c1 - c2 * x + c3 * x * w_h - c4 * x * w_h^2 + (c5 + c6 * x - c7 * x * w_h) / T,
+            c1 - c2 * x + c3 * x * w_h - c4 * x * w_h^2 +
+            (c5 + c6 * x - c7 * x * w_h) / T,
         ) * 100 # * 100 converts mbar --> Pa
     return p_sol
 end
@@ -202,7 +216,10 @@ end
 Returns the coefficients from Table B1 Appendix B in Chen et al 2022
 DOI: 10.1016/j.atmosres.2022.106171
 """
-function Chen2022_vel_coeffs_B1(coeffs::CMP.Chen2022VelTypeRain, ρₐ::FT) where {FT}
+function Chen2022_vel_coeffs_B1(
+    coeffs::CMP.Chen2022VelTypeRain,
+    ρₐ::FT,
+) where {FT}
     (; ρ0, a, a3_pow, b, b_ρ, c) = coeffs
     # Table B1
     q = exp(ρ0 * ρₐ)
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
index b4d9a40a4..f429aa844 100644
--- a/src/IceNucleation.jl
+++ b/src/IceNucleation.jl
@@ -200,7 +200,8 @@ function INP_concentration_frequency(
 
     μ = INP_concentration_mean(T)
 
-    return 1 / (sqrt(2 * FT(π)) * params.σ) * exp(-(log(INPC) - μ)^2 / (2 * params.σ^2))
+    return 1 / (sqrt(2 * FT(π)) * params.σ) *
+           exp(-(log(INPC) - μ)^2 / (2 * params.σ^2))
 end
 
 """
diff --git a/src/Microphysics0M.jl b/src/Microphysics0M.jl
index 4649ce5e6..64d9f1387 100644
--- a/src/Microphysics0M.jl
+++ b/src/Microphysics0M.jl
@@ -29,10 +29,15 @@ condensate specific humidity or supersaturation.
 The thresholds and the relaxation timescale are defined in
 ClimaParams.
 """
-remove_precipitation((; τ_precip, qc_0)::CMP.Parameters0M, q::TD.PhasePartition) =
-    -max(0, (q.liq + q.ice - qc_0)) / τ_precip
-
-remove_precipitation((; τ_precip, S_0)::CMP.Parameters0M, q::TD.PhasePartition, q_vap_sat) =
-    -max(0, (q.liq + q.ice - S_0 * q_vap_sat)) / τ_precip
+remove_precipitation(
+    (; τ_precip, qc_0)::CMP.Parameters0M,
+    q::TD.PhasePartition,
+) = -max(0, (q.liq + q.ice - qc_0)) / τ_precip
+
+remove_precipitation(
+    (; τ_precip, S_0)::CMP.Parameters0M,
+    q::TD.PhasePartition,
+    q_vap_sat,
+) = -max(0, (q.liq + q.ice - S_0 * q_vap_sat)) / τ_precip
 
 end #module Microphysics0M.jl
diff --git a/src/Microphysics1M.jl b/src/Microphysics1M.jl
index c14c59269..bc013d6d6 100644
--- a/src/Microphysics1M.jl
+++ b/src/Microphysics1M.jl
@@ -260,7 +260,8 @@ conv_q_liq_to_q_rai(
     q_liq::FT,
     smooth_transition::Bool = false,
 ) where {FT} =
-    smooth_transition ? CO.logistic_function_integral(q_liq, q_threshold, k) / τ :
+    smooth_transition ?
+    CO.logistic_function_integral(q_liq, q_threshold, k) / τ :
     max(0, q_liq - q_threshold) / τ
 
 """
@@ -280,7 +281,8 @@ conv_q_ice_to_q_sno_no_supersat(
     q_ice::FT,
     smooth_transition::Bool = false,
 ) where {FT} =
-    smooth_transition ? CO.logistic_function_integral(q_ice, q_threshold, k) / τ :
+    smooth_transition ?
+    CO.logistic_function_integral(q_ice, q_threshold, k) / τ :
     max(0, q_ice - q_threshold) / τ
 
 """
@@ -359,8 +361,8 @@ function accretion(
         E = Ec(cloud, precip, ce)
 
         accr_rate =
-            q_clo * E * n0 * a0 * v0 * χa * χv / λ * CO.Γ(ae + ve + Δa + Δv + FT(1)) /
-            (λ * r0)^(ae + ve + Δa + Δv)
+            q_clo * E * n0 * a0 * v0 * χa * χv / λ *
+            CO.Γ(ae + ve + Δa + Δv + FT(1)) / (λ * r0)^(ae + ve + Δa + Δv)
     end
     return accr_rate
 end
@@ -463,8 +465,10 @@ function accretion_snow_rain(
         accr_rate =
             FT(π) / ρ * n0_i * n0_j * m0_j * χm_j * E_ij * abs(v_ti - v_tj) /
             r0_j^(me_j + Δm_j) * (
-                FT(2) * CO.Γ(me_j + Δm_j + FT(1)) / λ_i^FT(3) / λ_j^(me_j + Δm_j + FT(1)) +
-                FT(2) * CO.Γ(me_j + Δm_j + FT(2)) / λ_i^FT(2) / λ_j^(me_j + Δm_j + FT(2)) +
+                FT(2) * CO.Γ(me_j + Δm_j + FT(1)) / λ_i^FT(3) /
+                λ_j^(me_j + Δm_j + FT(1)) +
+                FT(2) * CO.Γ(me_j + Δm_j + FT(2)) / λ_i^FT(2) /
+                λ_j^(me_j + Δm_j + FT(2)) +
                 CO.Γ(me_j + Δm_j + FT(3)) / λ_i / λ_j^(me_j + Δm_j + FT(3))
             )
     end
@@ -517,7 +521,8 @@ function evaporation_sublimation(
         evap_subl_rate =
             4 * FT(π) * n0 / ρ * S * G / λ^FT(2) * (
                 a_vent +
-                b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) *
+                b_vent * (ν_air / D_vapor)^FT(1 / 3) /
+                (r0 * λ)^((ve + Δv) / FT(2)) *
                 (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) *
                 CO.Γ((ve + Δv + FT(5)) / FT(2))
             )
@@ -555,7 +560,8 @@ function evaporation_sublimation(
         evap_subl_rate =
             4 * FT(π) * n0 / ρ * S * G / λ^FT(2) * (
                 a_vent +
-                b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) *
+                b_vent * (ν_air / D_vapor)^FT(1 / 3) /
+                (r0 * λ)^((ve + Δv) / FT(2)) *
                 (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) *
                 CO.Γ((ve + Δv + FT(5)) / FT(2))
             )
@@ -605,7 +611,8 @@ function snow_melt(
         snow_melt_rate =
             4 * FT(π) * n0 / ρ * K_therm / L * (T - T_freeze) / λ^FT(2) * (
                 a_vent +
-                b_vent * (ν_air / D_vapor)^FT(1 / 3) / (r0 * λ)^((ve + Δv) / FT(2)) *
+                b_vent * (ν_air / D_vapor)^FT(1 / 3) /
+                (r0 * λ)^((ve + Δv) / FT(2)) *
                 (FT(2) * v0 * χv / ν_air / λ)^FT(1 / 2) *
                 CO.Γ((ve + Δv + FT(5)) / FT(2))
             )
diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl
index 7aa237082..c0d07dfe0 100644
--- a/src/Microphysics2M.jl
+++ b/src/Microphysics2M.jl
@@ -171,7 +171,10 @@ end
     particle size
 """
 function size_distribution(
-    pdf::Union{CMP.RainParticlePDF_SB2006{FT}, CMP.RainParticlePDF_SB2006_limited{FT}},
+    pdf::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
     D::FT,
     q::FT,
     ρₐ::FT,
@@ -213,7 +216,10 @@ end
     of the bounds from numerical solutions.
 """
 function get_size_distribution_bound(
-    pdf::Union{CMP.RainParticlePDF_SB2006{FT}, CMP.RainParticlePDF_SB2006_limited{FT}},
+    pdf::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
     q::FT,
     N::FT,
     ρₐ::FT,
@@ -233,11 +239,15 @@ function get_size_distribution_bound(
     cloud_λ = pdf_cloud_parameters(pdf, q, ρₐ, N).Ec^(1 / ψc) * 1e3 # converting to m
     cloud_problem(x) =
         tolerance -
-        exp(-exp(x)^ψc * cloud_λ^ψc) *
-        (1 + exp(x)^ψc * cloud_λ^ψc + 1 / 2 * exp(x)^(2 * ψc) * cloud_λ^(2 * ψc))
+        exp(-exp(x)^ψc * cloud_λ^ψc) * (
+            1 +
+            exp(x)^ψc * cloud_λ^ψc +
+            1 / 2 * exp(x)^(2 * ψc) * cloud_λ^(2 * ψc)
+        )
     guess =
         log(0.5) +
-        (log(0.00025) - log(0.5)) / (log(1e12) - log(1e2)) * (log(cloud_λ^3) - log(10^2))
+        (log(0.00025) - log(0.5)) / (log(1e12) - log(1e2)) *
+        (log(cloud_λ^3) - log(10^2))
     log_cloud_x =
         RS.find_zero(
             cloud_problem,
@@ -366,7 +376,8 @@ function liquid_self_collection(
 
     L_liq = ρ * q_liq
 
-    dN_liq_dt_sc = -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * L_liq^2 - dN_liq_dt_au
+    dN_liq_dt_sc =
+        -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * L_liq^2 - dN_liq_dt_au
 
     return dN_liq_dt_sc
 end
@@ -409,7 +420,10 @@ Returns the raindrops number density tendency due to collisions of raindrops
 that produce larger raindrops (self-collection) for `scheme == SB2006Type`
 """
 function rain_self_collection(
-    pdf::Union{CMP.RainParticlePDF_SB2006{FT}, CMP.RainParticlePDF_SB2006_limited{FT}},
+    pdf::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
     self::CMP.SelfColSB2006{FT},
     q_rai::FT,
     ρ::FT,
@@ -443,7 +457,10 @@ Returns the raindrops number density tendency due to breakup of raindrops
 that produce smaller raindrops for `scheme == SB2006Type`
 """
 function rain_breakup(
-    pdf::Union{CMP.RainParticlePDF_SB2006{FT}, CMP.RainParticlePDF_SB2006_limited{FT}},
+    pdf::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
     brek::CMP.BreakupSB2006{FT},
     q_rai::FT,
     ρ::FT,
@@ -459,7 +476,8 @@ function rain_breakup(
     xr = pdf_rain_parameters(pdf, q_rai, ρ, N_rai).xr
     Dr = (xr * 6 / FT(π) / ρw)^FT(1 / 3)
     ΔD = Dr - Deq
-    phi_br = (Dr < Dr_th) ? FT(-1) : ((ΔD <= 0) ? kbr * ΔD : 2 * (exp(κbr * ΔD) - 1))
+    phi_br =
+        (Dr < Dr_th) ? FT(-1) : ((ΔD <= 0) ? kbr * ΔD : 2 * (exp(κbr * ΔD) - 1))
     dN_rai_dt_br = -(phi_br + 1) * dN_rai_dt_sc
 
     return dN_rai_dt_br
@@ -537,7 +555,8 @@ function rain_terminal_velocity(
     # TODO: Input argument list needs to be redesigned
 
     λr = pdf_rain_parameters(pdf_r, q_rai, ρ, N_rai).λr
-    _pa0, _pb0, _pa1, _pb1 = _sb_rain_terminal_velocity_helper(pdf_r, λr, aR, bR, cR)
+    _pa0, _pb0, _pa1, _pb1 =
+        _sb_rain_terminal_velocity_helper(pdf_r, λr, aR, bR, cR)
 
     vt0 =
         N_rai < eps(FT) ? FT(0) :
@@ -652,10 +671,13 @@ function rain_evaporation(
 
         t_star = (FT(6) * x_star / xr)^FT(1 / 3)
         a_vent_0 = av * Γ_incl(FT(-1), t_star) / FT(6)^FT(-2 / 3)
-        b_vent_0 = bv * Γ_incl(-FT(0.5) + FT(1.5) * β, t_star) / FT(6)^FT(β / 2 - FT(0.5))
+        b_vent_0 =
+            bv * Γ_incl(-FT(0.5) + FT(1.5) * β, t_star) /
+            FT(6)^FT(β / 2 - FT(0.5))
 
         a_vent_1 = av * SF.gamma(FT(2)) / FT(6)^FT(1 / 3)
-        b_vent_1 = bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
+        b_vent_1 =
+            bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
 
         N_Re = α * xr^β * sqrt(ρ0 / ρ) * Dr / ν_air
         Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
@@ -666,7 +688,8 @@ function rain_evaporation(
 
         # When xr = 0 evap_rate_0 becomes NaN. We replace NaN with 0 which is the limit of
         # evap_rate_0 for xr -> 0.
-        evap_rate_0 = N_rai < eps(FT) || xr / x_star < eps(FT) ? FT(0) : evap_rate_0
+        evap_rate_0 =
+            N_rai < eps(FT) || xr / x_star < eps(FT) ? FT(0) : evap_rate_0
         evap_rate_1 = q_rai < eps(FT) ? FT(0) : evap_rate_1
     end
 
@@ -751,7 +774,9 @@ function conv_q_liq_to_q_rai(
         return FT(0)
     else
         # Mean volume radius in microns (assuming spherical cloud droplets)
-        r_vol = (FT(3) * (q_liq * ρ) / FT(4) / FT(π) / ρ_w / N_d)^FT(1 / 3) * FT(1e6)
+        r_vol =
+            (FT(3) * (q_liq * ρ) / FT(4) / FT(π) / ρ_w / N_d)^FT(1 / 3) *
+            FT(1e6)
 
         # Assumed size distribution: modified gamma distribution
         β_6 = ((r_vol + FT(3)) / r_vol)^FT(1 / 3)
diff --git a/src/Nucleation.jl b/src/Nucleation.jl
index fc24f78b8..140a9b25f 100644
--- a/src/Nucleation.jl
+++ b/src/Nucleation.jl
@@ -48,7 +48,13 @@ Calculates the rate of binary H2SO4-H2O and ternary H2SO4-H2O-NH3 nucleation for
 The particle formation rate is parameterized using data from the CLOUD experiment, through neutral and ion-induced channels.
 This is an implementation of Dunne et al 1016 doi:10.1126/science.aaf2649 Appendix 8-10
 """
-function h2so4_nucleation_rate(h2so4_conc, nh3_conc, negative_ion_conc, temp, params)
+function h2so4_nucleation_rate(
+    h2so4_conc,
+    nh3_conc,
+    negative_ion_conc,
+    temp,
+    params,
+)
     # Change units from 1/m³ to 1/cm³
     h2so4_conc *= 1e-6
     nh3_conc *= 1e-6
@@ -59,7 +65,8 @@ function h2so4_nucleation_rate(h2so4_conc, nh3_conc, negative_ion_conc, temp, pa
     # Calculate factors
     k(T, u, v, w) = exp(u - exp(v * (T / 1000 - w)))
     f_y(h2so4, nh3, p_t_y, p_A_y, a_y) =
-        (nh3 / ref_conc) / (a_y + (h2so4 / ref_conc)^p_t_y / (nh3 / ref_conc)^p_A_y)
+        (nh3 / ref_conc) /
+        (a_y + (h2so4 / ref_conc)^p_t_y / (nh3 / ref_conc)^p_A_y)
     k_b_n = k(temp, params.u_b_n, params.v_b_n, params.w_b_n)
     k_b_i = k(temp, params.u_b_i, params.v_b_i, params.w_b_i)
     k_t_n = k(temp, params.u_t_n, params.v_t_n, params.w_t_n)
@@ -75,7 +82,10 @@ function h2so4_nucleation_rate(h2so4_conc, nh3_conc, negative_ion_conc, temp, pa
         k_t_n * f_n * (h2so4_conc / ref_conc)^params.p_t_n +
         k_t_i * f_i * (h2so4_conc / ref_conc)^params.p_t_i * negative_ion_conc
     # Convert to 1/m³/s
-    return (; binary_rate = binary_rate * 1e6, ternary_rate = ternary_rate * 1e6)
+    return (;
+        binary_rate = binary_rate * 1e6,
+        ternary_rate = ternary_rate * 1e6,
+    )
 end
 
 """
@@ -115,17 +125,29 @@ function organic_nucleation_rate(
             params.Y_MTO3 * k_MTO3 * monoterpene_conc * O3_conc +
             params.Y_MTOH * k_MTOH * monoterpene_conc * OH_conc
         ) / condensation_sink
-    return organic_nucleation_rate_hom_prescribed(negative_ion_conc, HOM_conc, params)
+    return organic_nucleation_rate_hom_prescribed(
+        negative_ion_conc,
+        HOM_conc,
+        params,
+    )
 end
 
-function organic_nucleation_rate_hom_prescribed(negative_ion_conc, HOM_conc, params)
+function organic_nucleation_rate_hom_prescribed(
+    negative_ion_conc,
+    HOM_conc,
+    params,
+)
     # HOM reference concentration: 1e7/cm³
     ref_conc = 1e7
     rate =
         params.a_1 *
-        (HOM_conc / ref_conc)^(params.a_2 + params.a_5 / (HOM_conc / ref_conc)) +
+        (
+            HOM_conc / ref_conc
+        )^(params.a_2 + params.a_5 / (HOM_conc / ref_conc)) +
         params.a_3 *
-        (HOM_conc / ref_conc)^(params.a_4 + params.a_5 / (HOM_conc / ref_conc)) *
+        (
+            HOM_conc / ref_conc
+        )^(params.a_4 + params.a_5 / (HOM_conc / ref_conc)) *
         negative_ion_conc
     # Convert from (1/cm³/s) to (1/m³/s)
     return rate * 1e6
@@ -163,7 +185,11 @@ function organic_and_h2so4_nucleation_rate(
     )
 end
 
-function organic_and_h2so4_nucleation_rate_bioOxOrg_prescribed(h2so4_conc, bioOxOrg, params)
+function organic_and_h2so4_nucleation_rate_bioOxOrg_prescribed(
+    h2so4_conc,
+    bioOxOrg,
+    params,
+)
     # Convert bioOxOrg and k_H2SO4org from 1/m³ to 1/cm³
     k_H2SO4org = 1e-6 * params.k_H2SO4org
     bioOxOrg *= 1e-6
diff --git a/src/P3_helpers.jl b/src/P3_helpers.jl
index 6794b9a65..af18ae512 100644
--- a/src/P3_helpers.jl
+++ b/src/P3_helpers.jl
@@ -71,7 +71,11 @@ function Γ_lower(a::FT, z::FT) where {FT <: Real}
 
         return exp(-z) *
                z^a *
-               (1 / a + c₁(a) * z / a / (a + 1) + (c₁(a) * z)^2 / a / (a + 1) / (a + 2)) *
+               (
+                   1 / a +
+                   c₁(a) * z / a / (a + 1) +
+                   (c₁(a) * z)^2 / a / (a + 1) / (a + 2)
+               ) *
                (1 - W(a, z)) + CO.Γ(a) * W(a, z) * (1 - c₄(a)^(-z))
     end
 end
@@ -102,7 +106,9 @@ function ∫_Γ(x₀::FT, x_end::FT, c1::FT, c2::FT, c3::FT) where {FT}
     elseif x₀ == 0
         return c1 * c3^(-c2 - 1) * (Γ_lower(1 + c2, x_end * c3))
     else
-        return c1 * c3^(-c2 - 1) * (Γ_upper(1 + c2, x₀ * c3) - Γ_upper(1 + c2, x_end * c3))
+        return c1 *
+               c3^(-c2 - 1) *
+               (Γ_upper(1 + c2, x₀ * c3) - Γ_upper(1 + c2, x_end * c3))
     end
 end
 
diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl
index 04606f26c..b4bdb66ec 100644
--- a/src/P3_particle_properties.jl
+++ b/src/P3_particle_properties.jl
@@ -62,7 +62,8 @@ end
 Returns the density of total (deposition + rime) ice mass for graupel, in kg/m3
 Eq. 16 in Morrison and Milbrandt (2015).
 """
-ρ_g_helper(ρ_r::FT, F_rim::FT, ρ_d::FT) where {FT} = F_rim * ρ_r + (1 - F_rim) * ρ_d
+ρ_g_helper(ρ_r::FT, F_rim::FT, ρ_d::FT) where {FT} =
+    F_rim * ρ_r + (1 - F_rim) * ρ_d
 
 """
     ρ_d_helper(p3, D_cr, D_gr)
@@ -77,7 +78,8 @@ Eq. 17 in Morrison and Milbrandt (2015).
 function ρ_d_helper(p3::PSP3{FT}, D_cr::FT, D_gr::FT) where {FT}
     α_va = α_va_si(p3)
     β_m2 = p3.β_va - 2
-    return 6 * α_va * (D_cr^β_m2 - D_gr^β_m2) / FT(π) / β_m2 / max(D_cr - D_gr, eps(FT))
+    return 6 * α_va * (D_cr^β_m2 - D_gr^β_m2) / FT(π) / β_m2 /
+           max(D_cr - D_gr, eps(FT))
 end
 
 """
@@ -121,7 +123,12 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_rim::FT) where {FT}
             ).root
         ρ_g = ρ_g_helper(ρ_r, F_rim, ρ_d)
 
-        return (; D_cr = D_cr_helper(p3, F_rim, ρ_g), D_gr = D_gr_helper(p3, ρ_g), ρ_g, ρ_d)
+        return (;
+            D_cr = D_cr_helper(p3, F_rim, ρ_g),
+            D_gr = D_gr_helper(p3, ρ_g),
+            ρ_g,
+            ρ_d,
+        )
     end
 end
 
@@ -262,7 +269,8 @@ A_s(D::FT) where {FT} = FT(π) / 4 * D^2
 # for nonspherical particles
 A_ns(p3::PSP3, D::FT) where {FT} = p3.γ * D^p3.σ
 # partially rimed ice
-A_r(p3::PSP3, F_rim::FT, D::FT) where {FT} = F_rim * A_s(D) + (1 - F_rim) * A_ns(p3, D)
+A_r(p3::PSP3, F_rim::FT, D::FT) where {FT} =
+    F_rim * A_s(D) + (1 - F_rim) * A_ns(p3, D)
 
 """
     p3_area(p3, D, F_rim, F_liq, th)
diff --git a/src/P3_processes.jl b/src/P3_processes.jl
index dd22ede00..515960a52 100644
--- a/src/P3_processes.jl
+++ b/src/P3_processes.jl
@@ -86,7 +86,8 @@ function ice_melt(
 
     # Get the P3 diameter distribution...
     th = thresholds(p3, ρ_rim, F_rim)
-    (λ, N_0) = distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_)
+    (λ, N_0) =
+        distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_)
     N(D) = N′ice(p3, D, λ, N_0)
     # ... and D_max for the integral
     bound = get_ice_bound(p3, λ, N_ice, FT(1e-6))
diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl
index e16f164a5..8da23715d 100644
--- a/src/P3_size_distribution.jl
+++ b/src/P3_size_distribution.jl
@@ -72,7 +72,8 @@ function get_ice_bound(p3::PSP3, λ::FT, N::FT, rtol::FT) where {FT}
     μ = DSD_μ(p3, λ)
     N₀ = DSD_N₀(p3, N, λ)
 
-    ice_bound_problem(D_max_hat) = N - N₀ / λ^(μ + 1) * Γ_lower(μ + 1, D_max_hat)
+    ice_bound_problem(D_max_hat) =
+        N - N₀ / λ^(μ + 1) * Γ_lower(μ + 1, D_max_hat)
 
     D_guess_low = FT(1e-6)
     D_guess_high = FT(1e-1)
@@ -189,7 +190,14 @@ end
     - zero if solving for ice core shape parameters
 Returns the approximated shape parameter μ for a given q and N value
 """
-function DSD_μ_approx(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_rim::FT, F_liq::FT) where {FT}
+function DSD_μ_approx(
+    p3::PSP3,
+    L::FT,
+    N::FT,
+    ρ_r::FT,
+    F_rim::FT,
+    F_liq::FT,
+) where {FT}
     # Get thresholds for given F_rim, ρ_r
     th = thresholds(p3, ρ_r, F_rim)
 
@@ -249,7 +257,8 @@ function get_bounds(
     Ll = L_over_N_gamma(p3, F_rim, F_liq, log(left), μ, th)
     Lr = L_over_N_gamma(p3, F_rim, F_liq, log(right), μ, th)
 
-    guess = left * (goal / (Ll))^((log(right) - log(left)) / (log(Lr) - log(Ll)))
+    guess =
+        left * (goal / (Ll))^((log(right) - log(left)) / (log(Lr) - log(Ll)))
 
     max = log(guess * exp(radius))
     min = log(guess)
@@ -342,7 +351,8 @@ function D_m(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_rim::FT, F_liq::FT) where {FT}
         n_nl += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ)
         n_nl += ∫_Γ(D_th, th.D_gr, α_va * N_0, μ + p3.β_va + 1, λ)
         n_nl += ∫_Γ(th.D_gr, th.D_cr, π / 6 * th.ρ_g * N_0, μ + 4, λ)
-        n_nl += ∫_Γ(th.D_cr, FT(Inf), α_va / (1 - F_rim) * N_0, μ + p3.β_va + 1, λ)
+        n_nl +=
+            ∫_Γ(th.D_cr, FT(Inf), α_va / (1 - F_rim) * N_0, μ + p3.β_va + 1, λ)
         n_l = ∫_Γ(FT(0), FT(Inf), N_0 * p3.ρ_l * (FT(π) / 6), μ + 4, λ)
     end
 
diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl
index 07bbe26d7..d089edc29 100644
--- a/src/P3_terminal_velocity.jl
+++ b/src/P3_terminal_velocity.jl
@@ -62,7 +62,15 @@ function p3_particle_terminal_velocity(
     th,
     use_aspect_ratio = true,
 ) where {FT}
-    v_i = ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th, use_aspect_ratio)
+    v_i = ice_particle_terminal_velocity(
+        p3,
+        D,
+        Chen2022,
+        ρₐ,
+        F_rim,
+        th,
+        use_aspect_ratio,
+    )
     v_r = CM2.rain_particle_terminal_velocity(D, Chen2022.rain, ρₐ)
     return p3_F_liq_average(F_liq, v_i, v_r)
 end
@@ -87,7 +95,10 @@ cloud or rain drop as a function of their sizes. It uses Chen 2022 velocity
 parameterization for ice and rain and assumes no sedimentation of cloud droplets.
 """
 function velocity_difference(
-    pdf_r::Union{CMP.RainParticlePDF_SB2006{FT}, CMP.RainParticlePDF_SB2006_limited{FT}},
+    pdf_r::Union{
+        CMP.RainParticlePDF_SB2006{FT},
+        CMP.RainParticlePDF_SB2006_limited{FT},
+    },
     Dₗ::FT,
     Dᵢ::FT,
     p3::PSP3,
diff --git a/src/PrecipitationSusceptibility.jl b/src/PrecipitationSusceptibility.jl
index 841bbfca5..fea1cebea 100644
--- a/src/PrecipitationSusceptibility.jl
+++ b/src/PrecipitationSusceptibility.jl
@@ -41,7 +41,9 @@ function precipitation_susceptibility_autoconversion(
     N_liq::FT,
 ) where {FT}
     grad = ForwardDiff.gradient(
-        x -> log(CM2.autoconversion(scheme.acnv, scheme.pdf_c, exp.(x)...).dq_rai_dt),
+        x -> log(
+            CM2.autoconversion(scheme.acnv, scheme.pdf_c, exp.(x)...).dq_rai_dt,
+        ),
         log.(abs.([q_liq, q_rai, ρ, N_liq])),
     )
     return precip_susceptibility_rates(
diff --git a/src/parameters/AerosolDesertDust.jl b/src/parameters/AerosolDesertDust.jl
index b3f954a4f..c413bb8f5 100644
--- a/src/parameters/AerosolDesertDust.jl
+++ b/src/parameters/AerosolDesertDust.jl
@@ -25,7 +25,8 @@ Base.@kwdef struct DesertDust{FT} <: AerosolType{FT}
     ABIFM_c::FT
 end
 
-DesertDust(::Type{FT}) where {FT <: AbstractFloat} = DesertDust(CP.create_toml_dict(FT))
+DesertDust(::Type{FT}) where {FT <: AbstractFloat} =
+    DesertDust(CP.create_toml_dict(FT))
 
 function DesertDust(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/AerosolFeldspar.jl b/src/parameters/AerosolFeldspar.jl
index eb7c939eb..6d072d9cd 100644
--- a/src/parameters/AerosolFeldspar.jl
+++ b/src/parameters/AerosolFeldspar.jl
@@ -16,7 +16,8 @@ Base.@kwdef struct Feldspar{FT} <: AerosolType{FT}
     deposition_c::FT
 end
 
-Feldspar(::Type{FT}) where {FT <: AbstractFloat} = Feldspar(CP.create_toml_dict(FT))
+Feldspar(::Type{FT}) where {FT <: AbstractFloat} =
+    Feldspar(CP.create_toml_dict(FT))
 
 function Feldspar(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/AerosolFerrihydrite.jl b/src/parameters/AerosolFerrihydrite.jl
index e9fa6709d..131fdc497 100644
--- a/src/parameters/AerosolFerrihydrite.jl
+++ b/src/parameters/AerosolFerrihydrite.jl
@@ -16,7 +16,8 @@ Base.@kwdef struct Ferrihydrite{FT} <: AerosolType{FT}
     deposition_c::FT
 end
 
-Ferrihydrite(::Type{FT}) where {FT <: AbstractFloat} = Ferrihydrite(CP.create_toml_dict(FT))
+Ferrihydrite(::Type{FT}) where {FT <: AbstractFloat} =
+    Ferrihydrite(CP.create_toml_dict(FT))
 
 function Ferrihydrite(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/AerosolKaolinite.jl b/src/parameters/AerosolKaolinite.jl
index e623c6e53..0b3a66b47 100644
--- a/src/parameters/AerosolKaolinite.jl
+++ b/src/parameters/AerosolKaolinite.jl
@@ -21,7 +21,8 @@ Base.@kwdef struct Kaolinite{FT} <: AerosolType{FT}
     ABIFM_c::FT
 end
 
-Kaolinite(::Type{FT}) where {FT <: AbstractFloat} = Kaolinite(CP.create_toml_dict(FT))
+Kaolinite(::Type{FT}) where {FT <: AbstractFloat} =
+    Kaolinite(CP.create_toml_dict(FT))
 
 function Kaolinite(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/AerosolModalNucleation.jl b/src/parameters/AerosolModalNucleation.jl
index 2e24488c3..564955f3f 100644
--- a/src/parameters/AerosolModalNucleation.jl
+++ b/src/parameters/AerosolModalNucleation.jl
@@ -1,4 +1,5 @@
-export H2S04NucleationParameters, OrganicNucleationParameters, MixedNucleationParameters
+export H2S04NucleationParameters,
+    OrganicNucleationParameters, MixedNucleationParameters
 
 """
     H2S04NucleationParameters{FT}
@@ -129,7 +130,8 @@ MixedNucleationParameters(::Type{FT}) where {FT <: AbstractFloat} =
 
 function MixedNucleationParameters(td::CP.AbstractTOMLDict)
     name_map = (;
-        :mam3_nucleation_k_H2SO4_mixed_organic_sulfuric_acid_factor => :k_H2SO4org,
+        :mam3_nucleation_k_H2SO4_mixed_organic_sulfuric_acid_factor =>
+            :k_H2SO4org,
         :mam3_nucleation_k_MTOH_organic_factor => :k_MTOH,
         :mam3_nucleation_exp_MTOH_organic_factor => :exp_MTOH,
     )
diff --git a/src/parameters/AerosolSeasalt.jl b/src/parameters/AerosolSeasalt.jl
index 152f2401a..bf92ae02e 100644
--- a/src/parameters/AerosolSeasalt.jl
+++ b/src/parameters/AerosolSeasalt.jl
@@ -23,7 +23,8 @@ Base.@kwdef struct Seasalt{FT} <: AerosolType{FT}
     κ::FT
 end
 
-Seasalt(::Type{FT}) where {FT <: AbstractFloat} = Seasalt(CP.create_toml_dict(FT))
+Seasalt(::Type{FT}) where {FT <: AbstractFloat} =
+    Seasalt(CP.create_toml_dict(FT))
 
 function Seasalt(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/AerosolSulfate.jl b/src/parameters/AerosolSulfate.jl
index 728008086..3b4375f2f 100644
--- a/src/parameters/AerosolSulfate.jl
+++ b/src/parameters/AerosolSulfate.jl
@@ -23,7 +23,8 @@ Base.@kwdef struct Sulfate{FT} <: AerosolType{FT}
     κ::FT
 end
 
-Sulfate(::Type{FT}) where {FT <: AbstractFloat} = Sulfate(CP.create_toml_dict(FT))
+Sulfate(::Type{FT}) where {FT <: AbstractFloat} =
+    Sulfate(CP.create_toml_dict(FT))
 
 function Sulfate(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl
index f1aae21a7..00052bac4 100644
--- a/src/parameters/IceNucleation.jl
+++ b/src/parameters/IceNucleation.jl
@@ -18,8 +18,10 @@ Base.@kwdef struct Mohler2006{FT} <: ParametersType{FT}
 end
 
 function Mohler2006(td::CP.AbstractTOMLDict)
-    name_map =
-        (; :Mohler2006_maximum_allowed_Si => :Sᵢ_max, :Mohler2006_threshold_T => :T_thr)
+    name_map = (;
+        :Mohler2006_maximum_allowed_Si => :Sᵢ_max,
+        :Mohler2006_threshold_T => :T_thr,
+    )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     FT = CP.float_type(td)
     return Mohler2006{FT}(; parameters...)
@@ -132,7 +134,11 @@ function IceNucleationParameters(toml_dict::CP.AbstractTOMLDict)
     HOM = typeof(homogeneous)
     P3_type = typeof(p3)
     FT = CP.float_type(toml_dict)
-    return IceNucleationParameters{FT, DEP, HOM, P3_type}(deposition, homogeneous, p3)
+    return IceNucleationParameters{FT, DEP, HOM, P3_type}(
+        deposition,
+        homogeneous,
+        p3,
+    )
 end
 
 
diff --git a/src/parameters/Microphysics0M.jl b/src/parameters/Microphysics0M.jl
index e9dc3504c..e704a095b 100644
--- a/src/parameters/Microphysics0M.jl
+++ b/src/parameters/Microphysics0M.jl
@@ -17,7 +17,8 @@ Base.@kwdef struct Parameters0M{FT} <: ParametersType{FT}
     S_0::FT
 end
 
-Parameters0M(::Type{FT}) where {FT <: AbstractFloat} = Parameters0M(CP.create_toml_dict(FT))
+Parameters0M(::Type{FT}) where {FT <: AbstractFloat} =
+    Parameters0M(CP.create_toml_dict(FT))
 
 function Parameters0M(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/Microphysics1M.jl b/src/parameters/Microphysics1M.jl
index f80cecd3a..8c5f9b5e4 100644
--- a/src/parameters/Microphysics1M.jl
+++ b/src/parameters/Microphysics1M.jl
@@ -137,7 +137,8 @@ Base.@kwdef struct CloudLiquid{FT} <: CloudCondensateType{FT}
     r_eff::FT
 end
 
-CloudLiquid(::Type{FT}) where {FT <: AbstractFloat} = CloudLiquid(CP.create_toml_dict(FT))
+CloudLiquid(::Type{FT}) where {FT <: AbstractFloat} =
+    CloudLiquid(CP.create_toml_dict(FT))
 
 function CloudLiquid(toml_dict::CP.AbstractTOMLDict)
     name_map = (;
@@ -145,7 +146,8 @@ function CloudLiquid(toml_dict::CP.AbstractTOMLDict)
         :density_liquid_water => :ρw,
         :liquid_cloud_effective_radius => :r_eff,
     )
-    parameters = CP.get_parameter_values(toml_dict, name_map, "CloudMicrophysics")
+    parameters =
+        CP.get_parameter_values(toml_dict, name_map, "CloudMicrophysics")
     FT = CP.float_type(toml_dict)
     return CloudLiquid{FT}(; parameters...)
 end
@@ -175,7 +177,8 @@ struct CloudIce{FT, PD, MS} <: CloudCondensateType{FT}
     r_eff::FT
 end
 
-CloudIce(::Type{FT}) where {FT <: AbstractFloat} = CloudIce(CP.create_toml_dict(FT))
+CloudIce(::Type{FT}) where {FT <: AbstractFloat} =
+    CloudIce(CP.create_toml_dict(FT))
 
 function CloudIce(toml_dict::CP.AbstractTOMLDict = CP.create_toml_dict(FT))
     name_map = (;
@@ -193,7 +196,15 @@ function CloudIce(toml_dict::CP.AbstractTOMLDict = CP.create_toml_dict(FT))
     FT = CP.float_type(toml_dict)
     P = typeof(pdf)
     M = typeof(mass)
-    return CloudIce{FT, P, M}(pdf, mass, p.r0, p.r_ice_snow, p.τ_relax, p.ρᵢ, p.r_eff)
+    return CloudIce{FT, P, M}(
+        pdf,
+        mass,
+        p.r0,
+        p.r_ice_snow,
+        p.τ_relax,
+        p.ρᵢ,
+        p.r_eff,
+    )
 end
 
 function ParticleMass(::Type{CloudIce}, td::CP.AbstractTOMLDict)
@@ -243,7 +254,8 @@ function Rain(toml_dict::CP.AbstractTOMLDict)
         :rain_mass_size_relation_coefficient_me => :me,
         :rain_cross_section_size_relation_coefficient_ae => :ae,
         :rain_autoconversion_timescale => :τ,
-        :cloud_liquid_water_specific_humidity_autoconversion_threshold => :q_threshold,
+        :cloud_liquid_water_specific_humidity_autoconversion_threshold =>
+            :q_threshold,
         :threshold_smooth_transition_steepness => :k,
         :rain_ventillation_coefficient_a => :a,
         :rain_ventillation_coefficient_b => :b,
@@ -328,7 +340,8 @@ function Snow(toml_dict::CP.AbstractTOMLDict)
         :cloud_ice_crystals_length_scale => :r0,
         :snow_apparent_density => :ρᵢ,
         :snow_autoconversion_timescale => :τ,
-        :cloud_ice_specific_humidity_autoconversion_threshold => :q_threshold,
+        :cloud_ice_specific_humidity_autoconversion_threshold =>
+            :q_threshold,
         :threshold_smooth_transition_steepness => :k,
         :temperature_water_freeze => :T_freeze,
         :snow_flake_size_distribution_coefficient_mu => :μ,
@@ -415,7 +428,8 @@ Base.@kwdef struct CollisionEff{FT} <: ParametersType{FT}
     e_rai_sno::FT
 end
 
-CollisionEff(::Type{FT}) where {FT <: AbstractFloat} = CollisionEff(CP.create_toml_dict(FT))
+CollisionEff(::Type{FT}) where {FT <: AbstractFloat} =
+    CollisionEff(CP.create_toml_dict(FT))
 
 function CollisionEff(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/Microphysics2M.jl b/src/parameters/Microphysics2M.jl
index f80051638..61d007836 100644
--- a/src/parameters/Microphysics2M.jl
+++ b/src/parameters/Microphysics2M.jl
@@ -144,8 +144,11 @@ struct AccrB1994{FT} <: ParametersType{FT}
 end
 
 function AccrB1994(toml_dict::CP.AbstractTOMLDict)
-    (; B1994_accretion_coeff_A) =
-        CP.get_parameter_values(toml_dict, "B1994_accretion_coeff_A", "CloudMicrophysics")
+    (; B1994_accretion_coeff_A) = CP.get_parameter_values(
+        toml_dict,
+        "B1994_accretion_coeff_A",
+        "CloudMicrophysics",
+    )
     FT = CP.float_type(toml_dict)
     return AccrB1994{FT}(B1994_accretion_coeff_A)
 end
@@ -231,8 +234,11 @@ struct AccrTC1980{FT} <: ParametersType{FT}
 end
 
 function AccrTC1980(toml_dict::CP.AbstractTOMLDict)
-    (; TC1980_accretion_coeff_A) =
-        CP.get_parameter_values(toml_dict, "TC1980_accretion_coeff_A", "CloudMicrophysics")
+    (; TC1980_accretion_coeff_A) = CP.get_parameter_values(
+        toml_dict,
+        "TC1980_accretion_coeff_A",
+        "CloudMicrophysics",
+    )
     FT = CP.float_type(toml_dict)
     return AccrTC1980{FT}(TC1980_accretion_coeff_A)
 end
diff --git a/src/parameters/MicrophysicsP3.jl b/src/parameters/MicrophysicsP3.jl
index 3432ead96..8e79cc773 100644
--- a/src/parameters/MicrophysicsP3.jl
+++ b/src/parameters/MicrophysicsP3.jl
@@ -38,7 +38,8 @@ Base.@kwdef struct ParametersP3{FT} <: ParametersType{FT}
     vent_b::FT
 end
 
-ParametersP3(::Type{FT}) where {FT <: AbstractFloat} = ParametersP3(CP.create_toml_dict(FT))
+ParametersP3(::Type{FT}) where {FT <: AbstractFloat} =
+    ParametersP3(CP.create_toml_dict(FT))
 
 function ParametersP3(td::CP.AbstractTOMLDict)
     name_map = (;
diff --git a/src/parameters/TerminalVelocity.jl b/src/parameters/TerminalVelocity.jl
index b4fedb4b7..87f74b9b0 100644
--- a/src/parameters/TerminalVelocity.jl
+++ b/src/parameters/TerminalVelocity.jl
@@ -91,7 +91,8 @@ struct Blk1MVelType{FT, R, S} <: TerminalVelocityType{FT}
     snow::S
 end
 
-Blk1MVelType(::Type{FT}) where {FT <: AbstractFloat} = Blk1MVelType(CP.create_toml_dict(FT))
+Blk1MVelType(::Type{FT}) where {FT <: AbstractFloat} =
+    Blk1MVelType(CP.create_toml_dict(FT))
 
 function Blk1MVelType(toml_dict::CP.AbstractTOMLDict)
     rain = Blk1MVelTypeRain(toml_dict)
@@ -272,7 +273,12 @@ function Chen2022VelType(toml_dict::CP.AbstractTOMLDict)
     small_ice = Chen2022VelTypeSmallIce(toml_dict)
     large_ice = Chen2022VelTypeLargeIce(toml_dict)
     FT = CP.float_type(toml_dict)
-    return Chen2022VelType{FT, typeof(rain), typeof(small_ice), typeof(large_ice)}(
+    return Chen2022VelType{
+        FT,
+        typeof(rain),
+        typeof(small_ice),
+        typeof(large_ice),
+    }(
         rain,
         small_ice,
         large_ice,
diff --git a/test/aerosol_activation_calibration.jl b/test/aerosol_activation_calibration.jl
index 77e8a14d9..63c26f321 100644
--- a/test/aerosol_activation_calibration.jl
+++ b/test/aerosol_activation_calibration.jl
@@ -37,7 +37,8 @@ function calibrate_ARG(
         mkdir(joinpath(fpath, "data"))
         mv(fname, joinpath(fpath, "data", fname), force = true)
     end
-    X_train, Y_train, initial_data = read_aerosol_dataset(joinpath(fpath, "data", fname))
+    X_train, Y_train, initial_data =
+        read_aerosol_dataset(joinpath(fpath, "data", fname))
 
     rng = Random.MersenneTwister(1)
 
@@ -54,11 +55,13 @@ function calibrate_ARG(
     all_mean_error_metrics = []
 
     for i in 1:N_samples
-        sample_inds = StatsBase.sample(1:DF.nrow(X_train), sample_size, replace = false)
+        sample_inds =
+            StatsBase.sample(1:DF.nrow(X_train), sample_size, replace = false)
         X_sample = X_train[sample_inds, :]
         Y_sample = Y_train[sample_inds]
 
-        initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble)
+        initial_ensemble =
+            EKP.construct_initial_ensemble(rng, prior, N_ensemble)
         ensemble_kalman_process = EKP.EnsembleKalmanProcess(
             initial_ensemble,
             Y_sample,
@@ -71,13 +74,25 @@ function calibrate_ARG(
 
         for j in 1:N_iterations_per_sample
             params_cur = EKP.get_ϕ_final(prior, ensemble_kalman_process)
-            errs = calibration_error_metrics(X_sample, Y_sample, params_cur, aip, tps, FT)
+            errs = calibration_error_metrics(
+                X_sample,
+                Y_sample,
+                params_cur,
+                aip,
+                tps,
+                FT,
+            )
             push!(mean_error_metrics, StatsBase.mean(errs))
 
             pred_ens = hcat(
                 [
-                    calibrated_prediction(X_sample, params_cur[:, k], aip, tps, FT) for
-                    k in 1:N_ensemble
+                    calibrated_prediction(
+                        X_sample,
+                        params_cur[:, k],
+                        aip,
+                        tps,
+                        FT,
+                    ) for k in 1:N_ensemble
                 ]...,
             )
 
@@ -85,7 +100,14 @@ function calibrate_ARG(
         end
 
         params_final = EKP.get_ϕ_final(prior, ensemble_kalman_process)
-        errs = calibration_error_metrics(X_sample, Y_sample, params_final, aip, tps, FT)
+        errs = calibration_error_metrics(
+            X_sample,
+            Y_sample,
+            params_final,
+            aip,
+            tps,
+            FT,
+        )
         push!(mean_error_metrics, StatsBase.mean(errs))
 
         params_mean_final = EKP.get_ϕ_mean_final(prior, ensemble_kalman_process)
@@ -129,9 +151,11 @@ function test_emulator(FT; rtols = [1e-4, 1e-3, 0.26], N_samples_calib = 2)
     ap_calib = CMP.AerosolActivationParameters(calib_params)
 
     TT.@test AA.N_activated_per_mode(ap_calib, ad, aip, tps, T, p, w, q)[1] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[1] rtol = rtols[1]
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[1] rtol =
+        rtols[1]
     TT.@test AA.N_activated_per_mode(ap_calib, ad, aip, tps, T, p, w, q)[2] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[2] rtol = rtols[2]
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[2] rtol =
+        rtols[2]
     for n in 1:N_samples_calib
         TT.@test errs[n][end] < rtols[3]
     end
diff --git a/test/aerosol_activation_emulators.jl b/test/aerosol_activation_emulators.jl
index f785394f3..ca0332813 100644
--- a/test/aerosol_activation_emulators.jl
+++ b/test/aerosol_activation_emulators.jl
@@ -40,13 +40,20 @@ function MLJFlux.build(builder::NNBuilder, rng, n_in, n_out)
     else
         push!(
             layers,
-            Flux.Dense(n_in => builder.layer_sizes[1], Flux.sigmoid_fast, init = init),
+            Flux.Dense(
+                n_in => builder.layer_sizes[1],
+                Flux.sigmoid_fast,
+                init = init,
+            ),
         )
     end
     for i in 1:num_hidden_layers
         push!(layers, Flux.Dropout(builder.dropout[i]))
         if i == num_hidden_layers
-            push!(layers, Flux.Dense(builder.layer_sizes[i] => n_out, init = init))
+            push!(
+                layers,
+                Flux.Dense(builder.layer_sizes[i] => n_out, init = init),
+            )
         else
             push!(
                 layers,
@@ -80,8 +87,12 @@ function MLJ.fit(model::GPRegressor, verbosity, X, y)
                 Distributions.pdf(Distributions.Normal(0.0, 0.5), x) for
                 x in X.mode_1_ARG_act_frac
             ])
-            inds =
-                StatsBase.sample(1:DF.nrow(X), weights, model.sample_size, replace = false)
+            inds = StatsBase.sample(
+                1:DF.nrow(X),
+                weights,
+                model.sample_size,
+                replace = false,
+            )
             inds_inducing = StatsBase.sample(
                 1:DF.nrow(X),
                 weights,
@@ -89,9 +100,16 @@ function MLJ.fit(model::GPRegressor, verbosity, X, y)
                 replace = false,
             )
         else
-            inds = StatsBase.sample(1:DF.nrow(X), model.sample_size, replace = false)
-            inds_inducing =
-                StatsBase.sample(1:DF.nrow(X), model.sample_size_inducing, replace = false)
+            inds = StatsBase.sample(
+                1:DF.nrow(X),
+                model.sample_size,
+                replace = false,
+            )
+            inds_inducing = StatsBase.sample(
+                1:DF.nrow(X),
+                model.sample_size_inducing,
+                replace = false,
+            )
         end
         if model.use_DTC
             gp = GaussianProcesses.DTC(
@@ -127,7 +145,13 @@ function MLJ.predict(::GPRegressor, fitresult, Xnew)
         hcat,
         [GaussianProcesses.predict_f(gp, Matrix(Xnew)')[2] for gp in fitresult],
     )
-    return (sum(means ./ variances, dims = 2) ./ sum(1.0 ./ variances, dims = 2))[:, 1]
+    return (sum(
+        means ./ variances,
+        dims = 2,
+    ) ./ sum(1.0 ./ variances, dims = 2))[
+        :,
+        1,
+    ]
 end
 
 # Load aerosol data reading and preprocessing functions
@@ -176,7 +200,8 @@ function get_2modal_model_FT32(;
 
         # Define the training pipeline
         if with_ARG_act_frac
-            _preprocess_aerosol_data = preprocess_aerosol_data_with_ARG_act_frac_FT32
+            _preprocess_aerosol_data =
+                preprocess_aerosol_data_with_ARG_act_frac_FT32
         else
             _preprocess_aerosol_data = preprocess_aerosol_data
         end
@@ -186,8 +211,15 @@ function get_2modal_model_FT32(;
                 Standardizer() |>
                 MLJ.TransformedTargetModel(
                     NeuralNetworkRegressor(
-                        builder = NNBuilder([250, 50, 5], [FT(0.3), FT(0), FT(0)]),
-                        optimiser = Flux.Optimise.Adam(0.001, (0.9, 0.999), 1.0e-8),
+                        builder = NNBuilder(
+                            [250, 50, 5],
+                            [FT(0.3), FT(0), FT(0)],
+                        ),
+                        optimiser = Flux.Optimise.Adam(
+                            0.001,
+                            (0.9, 0.999),
+                            1.0e-8,
+                        ),
                         epochs = 2000,
                         loss = Flux.mse,
                         batch_size = 1000,
@@ -221,7 +253,12 @@ function get_2modal_model_FT32(;
                                 upper = 1,
                                 scale = :log,
                             ),
-                            range(EvoTreeRegressor(), :max_depth, lower = 3, upper = 15),
+                            range(
+                                EvoTreeRegressor(),
+                                :max_depth,
+                                lower = 3,
+                                upper = 15,
+                            ),
                         ],
                     ),
                     transformer = target_transform,
@@ -234,7 +271,9 @@ function get_2modal_model_FT32(;
         # Create the untrained ML model
         mach = MLJ.machine(pipeline, X_train, Y_train)
         # Train a new ML model
-        @info("Training a new ML model ($(ml_model)). This may take a minute or two.")
+        @info(
+            "Training a new ML model ($(ml_model)). This may take a minute or two."
+        )
         MLJ.fit!(mach, verbosity = 0)
         # Save the ML model to a binary file that can be re-used by CloudMicrophysics
         MLJ.save(joinpath(fpath, machine_name), mach)
@@ -284,9 +323,11 @@ function test_emulator(
     )
 
     TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q)[1] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[1] rtol = rtols[1]
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[1] rtol =
+        rtols[1]
     TT.@test AA.N_activated_per_mode(mach, ap, ad, aip, tps, T, p, w, q)[2] ≈
-             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[2] rtol = rtols[2]
+             AA.N_activated_per_mode(ap, ad, aip, tps, T, p, w, q)[2] rtol =
+        rtols[2]
     TT.@test AA.total_N_activated(mach, ap, ad, aip, tps, T, p, w, q) ≈
              AA.total_N_activated(ap, ad, aip, tps, T, p, w, q) rtol = rtols[2]
 
diff --git a/test/aerosol_activation_tests.jl b/test/aerosol_activation_tests.jl
index 75ad51072..a9fc9a559 100644
--- a/test/aerosol_activation_tests.jl
+++ b/test/aerosol_activation_tests.jl
@@ -17,7 +17,8 @@ function test_aerosol_activation(FT)
     #default ARG2000 parameter values
     ap_default = CMP.AerosolActivationParameters(FT)
     # calibrated values based on PySDM
-    override_file = joinpath(pkgdir(CM), "src", "parameters", "toml", "ARG2000.toml")
+    override_file =
+        joinpath(pkgdir(CM), "src", "parameters", "toml", "ARG2000.toml")
     toml_dict = CP.create_toml_dict(FT; override_file)
     ap_calibrated = CMP.AerosolActivationParameters(toml_dict)
 
@@ -133,11 +134,28 @@ function test_aerosol_activation(FT)
         for AM_t in (AM_3_B, AM_3_κ)
             for ap in (ap_default, ap_calibrated)
                 TT.@test all(AA.mean_hygroscopicity_parameter(ap, AM_t) .> 0.0)
-                TT.@test AA.max_supersaturation(ap, AM_t, aip, tps, T, p, w, q) > 0.0
-                TT.@test all(AA.N_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> 0.0)
-                TT.@test all(AA.M_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> 0.0)
-                TT.@test AA.total_N_activated(ap, AM_t, aip, tps, T, p, w, q) > 0.0
-                TT.@test AA.total_M_activated(ap, AM_t, aip, tps, T, p, w, q) > 0.0
+                TT.@test AA.max_supersaturation(
+                    ap,
+                    AM_t,
+                    aip,
+                    tps,
+                    T,
+                    p,
+                    w,
+                    q,
+                ) > 0.0
+                TT.@test all(
+                    AA.N_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .>
+                    0.0,
+                )
+                TT.@test all(
+                    AA.M_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .>
+                    0.0,
+                )
+                TT.@test AA.total_N_activated(ap, AM_t, aip, tps, T, p, w, q) >
+                         0.0
+                TT.@test AA.total_M_activated(ap, AM_t, aip, tps, T, p, w, q) >
+                         0.0
             end
         end
     end
@@ -237,9 +255,11 @@ function test_aerosol_activation(FT)
                 AD_κ = AM.AerosolDistribution((paper_mode_1_κ, paper_mode_2_κ))
 
                 N_act_frac_B[it] =
-                    AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / N_1_paper
+                    AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] /
+                    N_1_paper
                 N_act_frac_κ[it] =
-                    AA.N_activated_per_mode(ap, AD_κ, aip, tps, T, p, w, q)[1] / N_1_paper
+                    AA.N_activated_per_mode(ap, AD_κ, aip, tps, T, p, w, q)[1] /
+                    N_1_paper
                 it += 1
             end
 
diff --git a/test/aqua.jl b/test/aqua.jl
index f174dd5d9..3ab126849 100644
--- a/test/aqua.jl
+++ b/test/aqua.jl
@@ -22,7 +22,10 @@ using Aqua
         ambs = Aqua.detect_ambiguities(CloudMicrophysics; recursive = true)
         pkg_match(pkgname, pkdir::Nothing) = false
         pkg_match(pkgname, pkdir::AbstractString) = occursin(pkgname, pkdir)
-        filter!(x -> pkg_match("CloudMicrophysics", pkgdir(last(x).module)), ambs)
+        filter!(
+            x -> pkg_match("CloudMicrophysics", pkgdir(last(x).module)),
+            ambs,
+        )
 
         # If the number of ambiguities is less than the limit below,
         # then please lower the limit based on the new number of ambiguities.
diff --git a/test/cloud_diagnostics.jl b/test/cloud_diagnostics.jl
index fca8e8409..596bd622a 100644
--- a/test/cloud_diagnostics.jl
+++ b/test/cloud_diagnostics.jl
@@ -11,8 +11,13 @@ import CloudMicrophysics.CloudDiagnostics as CMD
 function test_cloud_diagnostics(FT)
 
     # Seifert and Beheng 2006 parameters
-    override_file =
-        joinpath(pkgdir(CM), "src", "parameters", "toml", "SB2006_limiters.toml")
+    override_file = joinpath(
+        pkgdir(CM),
+        "src",
+        "parameters",
+        "toml",
+        "SB2006_limiters.toml",
+    )
     toml_dict = CP.create_toml_dict(FT; override_file)
     SB2006 = CMP.SB2006(toml_dict)
     SB2006_no_limiters = CMP.SB2006(toml_dict, false)
@@ -29,11 +34,13 @@ function test_cloud_diagnostics(FT)
         ρ_air = FT(1)
         q_rai = FT(0.18e-3)
 
-        TT.@test CMD.radar_reflectivity_1M(rain, q_rai, ρ_air) ≈ FT(12.17) atol = 0.2
+        TT.@test CMD.radar_reflectivity_1M(rain, q_rai, ρ_air) ≈ FT(12.17) atol =
+            0.2
 
         q_rai = FT(0.89e-4)
 
-        TT.@test CMD.radar_reflectivity_1M(rain, q_rai, ρ_air) ≈ FT(6.68) atol = 0.2
+        TT.@test CMD.radar_reflectivity_1M(rain, q_rai, ρ_air) ≈ FT(6.68) atol =
+            0.2
 
     end
 
@@ -74,7 +81,14 @@ function test_cloud_diagnostics(FT)
         N_rai = FT(510859)
 
         #action
-        reff = CMD.effective_radius_Liu_Hallet_97(wtr, ρ_air, q_liq, N_liq, q_rai, N_rai)
+        reff = CMD.effective_radius_Liu_Hallet_97(
+            wtr,
+            ρ_air,
+            q_liq,
+            N_liq,
+            q_rai,
+            N_rai,
+        )
         #test
         TT.@test reff ≈ FT(2.66e-05) atol = FT(8e-6)
 
diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl
index e6fa9337b..adf3b0754 100644
--- a/test/common_functions_tests.jl
+++ b/test/common_functions_tests.jl
@@ -17,7 +17,11 @@ TT.@testset "logistic_function unit tests" begin
 
     TT.@test CO.logistic_function(1.0, 0.0, 2.0) == 1.0
     TT.@test CO.logistic_function(0.0, 0.0, 2.0) == 0.0
-    TT.@test_throws AssertionError("x_0 >= 0") CO.logistic_function(1.0, -1.0, 2.0)
+    TT.@test_throws AssertionError("x_0 >= 0") CO.logistic_function(
+        1.0,
+        -1.0,
+        2.0,
+    )
     TT.@test_throws AssertionError("k > 0") CO.logistic_function(1.0, 1.0, 0.0)
 
 end
@@ -31,8 +35,16 @@ TT.@testset "logistic_function_integral unit tests" begin
 
     TT.@test CO.logistic_function_integral(1.0, 0.0, 2.0) == 1.0
     TT.@test CO.logistic_function_integral(0.0, 0.0, 2.0) == 0.0
-    TT.@test_throws AssertionError("x_0 >= 0") CO.logistic_function_integral(1.0, -1.0, 2.0)
-    TT.@test_throws AssertionError("k > 0") CO.logistic_function_integral(1.0, 1.0, 0.0)
+    TT.@test_throws AssertionError("x_0 >= 0") CO.logistic_function_integral(
+        1.0,
+        -1.0,
+        2.0,
+    )
+    TT.@test_throws AssertionError("k > 0") CO.logistic_function_integral(
+        1.0,
+        1.0,
+        0.0,
+    )
 
 end
 
@@ -61,8 +73,15 @@ function test_H2SO4_soln_saturation_vapor_pressure(FT)
         )
 
         # p_sol should be higher at warmer temperatures
-        TT.@test CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph, T_warm) >
-                 CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph, T_cold)
+        TT.@test CO.H2SO4_soln_saturation_vapor_pressure(
+            H2SO4_prs,
+            x_sulph,
+            T_warm,
+        ) > CO.H2SO4_soln_saturation_vapor_pressure(
+            H2SO4_prs,
+            x_sulph,
+            T_cold,
+        )
     end
 end
 
@@ -141,27 +160,51 @@ function test_Chen_coefficients(FT)
         TT.@test all(
             isapprox.(
                 aiu,
-                [FT(286768.02047954104), FT(-1.6916433443360287e6), FT(9843.240767655458)],
+                [
+                    FT(286768.02047954104),
+                    FT(-1.6916433443360287e6),
+                    FT(9843.240767655458),
+                ],
+                rtol = tol,
+            ),
+        )
+        TT.@test all(
+            isapprox.(
+                bi,
+                [FT(2.249342), FT(2.249342), FT(1.098942)],
                 rtol = tol,
             ),
         )
-        TT.@test all(isapprox.(bi, [FT(2.249342), FT(2.249342), FT(1.098942)], rtol = tol))
-        TT.@test all(isapprox.(ciu, [FT(0), FT(184.325), FT(184.325)], rtol = tol))
+        TT.@test all(
+            isapprox.(ciu, [FT(0), FT(184.325), FT(184.325)], rtol = tol),
+        )
     end
 
     TT.@testset "Chen terminal velocity small ice (B2)" begin
         aiu, bi, ciu = CO.Chen2022_vel_coeffs_B2(Ch2022.small_ice, ρ, ice.ρᵢ)
 
-        TT.@test all(isapprox.(aiu, [312.9777159510928, -316.5335670126842], rtol = tol))
-        TT.@test all(isapprox.(bi, [0.7295470725655279, 0.7295470725655279], rtol = tol))
+        TT.@test all(
+            isapprox.(aiu, [312.9777159510928, -316.5335670126842], rtol = tol),
+        )
+        TT.@test all(
+            isapprox.(bi, [0.7295470725655279, 0.7295470725655279], rtol = tol),
+        )
         TT.@test all(isapprox.(ciu, [0.0, 4715.089121981011], rtol = tol))
     end
 
     TT.@testset "Chen terminal velocity large ice (B4)" begin
         aiu, bi, ciu = CO.Chen2022_vel_coeffs_B4(Ch2022.large_ice, ρ, snow.ρᵢ)
 
-        TT.@test all(isapprox.(aiu, [51.86069839334009, -1.394567234046072], rtol = tol))
-        TT.@test all(isapprox.(bi, [0.5655671081749194, 0.18155881980108224], rtol = tol))
+        TT.@test all(
+            isapprox.(aiu, [51.86069839334009, -1.394567234046072], rtol = tol),
+        )
+        TT.@test all(
+            isapprox.(
+                bi,
+                [0.5655671081749194, 0.18155881980108224],
+                rtol = tol,
+            ),
+        )
         TT.@test all(isapprox.(ciu, [0.0, 34.820462392120504], rtol = tol))
     end
 end
diff --git a/test/gpu_clima_core_test.jl b/test/gpu_clima_core_test.jl
index 64df5a066..f2debc2e7 100644
--- a/test/gpu_clima_core_test.jl
+++ b/test/gpu_clima_core_test.jl
@@ -76,7 +76,12 @@ function make_extruded_sphere(::Type{FT}) where {FT}
 
     # Define grid
     deep = false
-    grid = CC.Grids.ExtrudedFiniteDifferenceGrid(horz_grid, vert_grid, hypsography; deep)
+    grid = CC.Grids.ExtrudedFiniteDifferenceGrid(
+        horz_grid,
+        vert_grid,
+        hypsography;
+        deep,
+    )
 
     # Define 3D space
     center_extruded_space = CC.Spaces.CenterExtrudedFiniteDifferenceSpace(grid)
@@ -91,7 +96,12 @@ function set_sedimentation_precomputed_quantities(Y, p, t)
     (; w) = p
     (; params) = p
 
-    @. w = CMN.terminal_velocity(params.liquid, params.Ch2022.rain, Y.ρ, max(0, Y.ρq / Y.ρ))
+    @. w = CMN.terminal_velocity(
+        params.liquid,
+        params.Ch2022.rain,
+        Y.ρ,
+        max(0, Y.ρq / Y.ρ),
+    )
     return nothing
 end
 
diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl
index 2fb758238..182f8db82 100644
--- a/test/gpu_tests.jl
+++ b/test/gpu_tests.jl
@@ -77,7 +77,15 @@ const ArrayType = CuArray
             (ν[i],),
             (ρ[i],),
         )
-        mode_κ = AM.Mode_κ(r[i], stdev[i], N[i], (FT(1.0),), (FT(1.0),), (M[i],), (κ[i],))
+        mode_κ = AM.Mode_κ(
+            r[i],
+            stdev[i],
+            N[i],
+            (FT(1.0),),
+            (FT(1.0),),
+            (M[i],),
+            (κ[i],),
+        )
 
         arsl_dst_B = AM.AerosolDistribution((mode_B,))
         arsl_dst_κ = AM.AerosolDistribution((mode_κ,))
@@ -141,9 +149,11 @@ end
 
     @inbounds begin
         output[1, i] = CMN.terminal_velocity(liquid, Ch2022.rain, ρₐ[i], qₗ[i])
-        output[2, i] = CMN.terminal_velocity(ice, Ch2022.small_ice, ρₐ[i], qᵢ[i])
+        output[2, i] =
+            CMN.terminal_velocity(ice, Ch2022.small_ice, ρₐ[i], qᵢ[i])
         output[3, i] = CM1.terminal_velocity(rain, Ch2022.rain, ρₐ[i], qᵣ[i])
-        output[4, i] = CM1.terminal_velocity(snow, Ch2022.large_ice, ρₐ[i], qₛ[i])
+        output[4, i] =
+            CM1.terminal_velocity(snow, Ch2022.large_ice, ρₐ[i], qₛ[i])
     end
 end
 
@@ -204,12 +214,23 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[1, i] = CM1.accretion(liquid, rain, blk1mvel.rain, ce, ql[i], qr[i], ρ[i])
-        output[2, i] = CM1.accretion(ice, snow, blk1mvel.snow, ce, qi[i], qs[i], ρ[i])
-        output[3, i] = CM1.accretion(liquid, snow, blk1mvel.snow, ce, ql[i], qs[i], ρ[i])
-        output[4, i] = CM1.accretion(ice, rain, blk1mvel.rain, ce, qi[i], qr[i], ρ[i])
-        output[5, i] =
-            CM1.accretion_rain_sink(rain, ice, blk1mvel.rain, ce, qi[i], qr[i], ρ[i])
+        output[1, i] =
+            CM1.accretion(liquid, rain, blk1mvel.rain, ce, ql[i], qr[i], ρ[i])
+        output[2, i] =
+            CM1.accretion(ice, snow, blk1mvel.snow, ce, qi[i], qs[i], ρ[i])
+        output[3, i] =
+            CM1.accretion(liquid, snow, blk1mvel.snow, ce, ql[i], qs[i], ρ[i])
+        output[4, i] =
+            CM1.accretion(ice, rain, blk1mvel.rain, ce, qi[i], qr[i], ρ[i])
+        output[5, i] = CM1.accretion_rain_sink(
+            rain,
+            ice,
+            blk1mvel.rain,
+            ce,
+            qi[i],
+            qr[i],
+            ρ[i],
+        )
         output[6, i] = CM1.accretion_snow_rain(
             snow,
             rain,
@@ -247,7 +268,8 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[i] = CM1.snow_melt(snow, blk1mvel.snow, aps, tps, qs[i], ρ[i], T[i])
+        output[i] =
+            CM1.snow_melt(snow, blk1mvel.snow, aps, tps, qs[i], ρ[i], T[i])
     end
 end
 
@@ -354,18 +376,44 @@ end
                 ρ[i],
                 Nl[i],
             ).sc
-        output[6, i] = CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dq_liq_dt
-        output[7, i] = CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dN_liq_dt
-        output[8, i] = CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dq_rai_dt
-        output[9, i] = CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dN_rai_dt
-        output[10, i] = CM2.rain_self_collection_and_breakup(SB2006, qr[i], ρ[i], Nr[i]).sc
-        output[11, i] = CM2.rain_self_collection_and_breakup(SB2006, qr[i], ρ[i], Nr[i]).br
-        output[12, i] = CM2.rain_terminal_velocity(SB2006, SB2006Vel, qr[i], ρ[i], Nr[i])[1]
-        output[13, i] = CM2.rain_terminal_velocity(SB2006, SB2006Vel, qr[i], ρ[i], Nr[i])[2]
+        output[6, i] =
+            CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dq_liq_dt
+        output[7, i] =
+            CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dN_liq_dt
+        output[8, i] =
+            CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dq_rai_dt
+        output[9, i] =
+            CM2.accretion(SB2006, ql[i], qr[i], ρ[i], Nl[i]).dN_rai_dt
+        output[10, i] =
+            CM2.rain_self_collection_and_breakup(SB2006, qr[i], ρ[i], Nr[i]).sc
+        output[11, i] =
+            CM2.rain_self_collection_and_breakup(SB2006, qr[i], ρ[i], Nr[i]).br
+        output[12, i] =
+            CM2.rain_terminal_velocity(SB2006, SB2006Vel, qr[i], ρ[i], Nr[i])[1]
+        output[13, i] =
+            CM2.rain_terminal_velocity(SB2006, SB2006Vel, qr[i], ρ[i], Nr[i])[2]
         output[14, i] =
-            CM2.rain_evaporation(SB2006, aps, tps, q, qr[i], ρ[i], Nr[i], T[i]).evap_rate_0
+            CM2.rain_evaporation(
+                SB2006,
+                aps,
+                tps,
+                q,
+                qr[i],
+                ρ[i],
+                Nr[i],
+                T[i],
+            ).evap_rate_0
         output[15, i] =
-            CM2.rain_evaporation(SB2006, aps, tps, q, qr[i], ρ[i], Nr[i], T[i]).evap_rate_1
+            CM2.rain_evaporation(
+                SB2006,
+                aps,
+                tps,
+                q,
+                qr[i],
+                ρ[i],
+                Nr[i],
+                T[i],
+            ).evap_rate_1
     end
 end
 
@@ -477,7 +525,8 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[i] = CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph[i], T[i])
+        output[i] =
+            CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph[i], T[i])
     end
 end
 
@@ -496,7 +545,12 @@ end
     end
 end
 
-@kernel function Common_a_w_eT_kernel!(tps, output::AbstractArray{FT}, e, T) where {FT}
+@kernel function Common_a_w_eT_kernel!(
+    tps,
+    output::AbstractArray{FT},
+    e,
+    T,
+) where {FT}
 
     i = @index(Group, Linear)
 
@@ -505,7 +559,11 @@ end
     end
 end
 
-@kernel function Common_a_w_ice_kernel!(tps, output::AbstractArray{FT}, T) where {FT}
+@kernel function Common_a_w_ice_kernel!(
+    tps,
+    output::AbstractArray{FT},
+    T,
+) where {FT}
 
     i = @index(Group, Linear)
 
@@ -526,8 +584,14 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[1] = CMI_het.dust_activated_number_fraction(desert_dust, ip, S_i, T)
-        output[2] = CMI_het.dust_activated_number_fraction(arizona_test_dust, ip, S_i, T)
+        output[1] =
+            CMI_het.dust_activated_number_fraction(desert_dust, ip, S_i, T)
+        output[2] = CMI_het.dust_activated_number_fraction(
+            arizona_test_dust,
+            ip,
+            S_i,
+            T,
+        )
     end
 end
 
@@ -545,9 +609,22 @@ end
     i = @index(Group, Linear)
 
     @inbounds begin
-        output[1] = CMI_het.MohlerDepositionRate(desert_dust, ip, S_i, T, dSi_dt, N_aero)
-        output[2] =
-            CMI_het.MohlerDepositionRate(arizona_test_dust, ip, S_i, T, dSi_dt, N_aero)
+        output[1] = CMI_het.MohlerDepositionRate(
+            desert_dust,
+            ip,
+            S_i,
+            T,
+            dSi_dt,
+            N_aero,
+        )
+        output[2] = CMI_het.MohlerDepositionRate(
+            arizona_test_dust,
+            ip,
+            S_i,
+            T,
+            dSi_dt,
+            N_aero,
+        )
     end
 end
 
@@ -652,7 +729,12 @@ end
     end
 end
 
-@kernel function P3_scheme_kernel!(p3, output::AbstractArray{FT}, F_rim, ρ_r) where {FT}
+@kernel function P3_scheme_kernel!(
+    p3,
+    output::AbstractArray{FT},
+    F_rim,
+    ρ_r,
+) where {FT}
 
     i = @index(Group, Linear)
 
@@ -744,10 +826,16 @@ function test_gpu(FT)
         # test if all aerosol activation output is positive
         @test all(Array(output)[:, :] .>= FT(0))
         # test if higroscopicity parameter is the same for κ and B modes
-        @test all(isapprox(Array(output)[1, :], Array(output)[2, :], rtol = 0.3))
+        @test all(
+            isapprox(Array(output)[1, :], Array(output)[2, :], rtol = 0.3),
+        )
         # test if the number and mass activated are the same for κ and B modes
-        @test all(isapprox(Array(output)[3, :], Array(output)[4, :], rtol = 1e-5))
-        @test all(isapprox(Array(output)[5, :], Array(output)[6, :], rtol = 1e-5))
+        @test all(
+            isapprox(Array(output)[3, :], Array(output)[4, :], rtol = 1e-5),
+        )
+        @test all(
+            isapprox(Array(output)[5, :], Array(output)[6, :], rtol = 1e-5),
+        )
     end
 
     @testset "non-equilibrium microphysics kernels" begin
@@ -779,7 +867,20 @@ function test_gpu(FT)
         qₛ = ArrayType([FT(0.001)])
 
         kernel! = test_chen2022_terminal_velocity_kernel!(backend, work_groups)
-        kernel!(liquid, ice, rain, snow, Ch2022, output, ρ, qₗ, qᵢ, qᵣ, qₛ; ndrange)
+        kernel!(
+            liquid,
+            ice,
+            rain,
+            snow,
+            Ch2022,
+            output,
+            ρ,
+            qₗ,
+            qᵢ,
+            qᵣ,
+            qₛ;
+            ndrange,
+        )
 
         # test that terminal velocity is callable and returns a reasonable value
         @test Array(output)[1] ≈ FT(0.00689286343412659)
@@ -808,7 +909,8 @@ function test_gpu(FT)
         (; output, ndrange) = setup_output(dims, FT)
         ql = ArrayType([FT(1e-3), FT(5e-4)])
 
-        kernel! = test_1_moment_micro_smooth_transition_kernel!(backend, work_groups)
+        kernel! =
+            test_1_moment_micro_smooth_transition_kernel!(backend, work_groups)
         kernel!(rain.acnv1M, output, ql; ndrange)
 
         # Sanity checks for the GPU KernelAbstractions workflow
@@ -880,7 +982,18 @@ function test_gpu(FT)
         Nd = ArrayType([FT(1e8)])
 
         kernel! = test_2_moment_acnv_kernel!(backend, work_groups)
-        kernel!(KK2000, B1994, TC1980, LD2004, VarTSc, output, ql, ρ, Nd, ndrange = ndrange)
+        kernel!(
+            KK2000,
+            B1994,
+            TC1980,
+            LD2004,
+            VarTSc,
+            output,
+            ql,
+            ρ,
+            Nd,
+            ndrange = ndrange,
+        )
 
         @test Array(output)[1] ≈ FT(2e-6)
         @test Array(output)[2] ≈ FT(1.6963072465911614e-6)
@@ -946,8 +1059,16 @@ function test_gpu(FT)
                 @test isapprox(Array(output)[11], FT(14154.027), rtol = 1e-6)
                 @test isapprox(Array(output)[12], FT(0.9868878), rtol = 1e-6)
                 @test isapprox(Array(output)[13], FT(4.517734), rtol = 1e-6)
-                @test isapprox(Array(output)[14], FT(-243259.75126), rtol = 1e-6)
-                @test isapprox(Array(output)[15], FT(-0.0034601581), rtol = 1e-6)
+                @test isapprox(
+                    Array(output)[14],
+                    FT(-243259.75126),
+                    rtol = 1e-6,
+                )
+                @test isapprox(
+                    Array(output)[15],
+                    FT(-0.0034601581),
+                    rtol = 1e-6,
+                )
             end
             if SB == SB2006_no_limiters
                 @test isapprox(Array(output)[10], FT(-40447.855), rtol = 1e-6)
@@ -955,7 +1076,11 @@ function test_gpu(FT)
                 @test isapprox(Array(output)[12], FT(2.6429e-3), rtol = 1e-4)
                 @test isapprox(Array(output)[13], FT(0.1149338), rtol = 1e-6)
                 @test isapprox(Array(output)[14], FT(-52903.817), rtol = 1e-6)
-                @test isapprox(Array(output)[15], FT(-9.3601206e-5), rtol = 1e-6)
+                @test isapprox(
+                    Array(output)[15],
+                    FT(-9.3601206e-5),
+                    rtol = 1e-6,
+                )
             end
         end
     end
@@ -967,7 +1092,10 @@ function test_gpu(FT)
         T = ArrayType([FT(230)])
         x_sulph = ArrayType([FT(0.1)])
 
-        kernel! = Common_H2SO4_soln_saturation_vapor_pressure_kernel!(backend, work_groups)
+        kernel! = Common_H2SO4_soln_saturation_vapor_pressure_kernel!(
+            backend,
+            work_groups,
+        )
         kernel!(H2SO4_prs, output, x_sulph, T; ndrange)
 
         # test H2SO4_soln_saturation_vapor_pressure is callable and returns a reasonable value
@@ -1016,8 +1144,19 @@ function test_gpu(FT)
         T = FT(240)
         S_i = FT(1.2)
 
-        kernel! = IceNucleation_dust_activated_number_fraction_kernel!(backend, work_groups)
-        kernel!(output, desert_dust, arizona_test_dust, ip.deposition, S_i, T; ndrange)
+        kernel! = IceNucleation_dust_activated_number_fraction_kernel!(
+            backend,
+            work_groups,
+        )
+        kernel!(
+            output,
+            desert_dust,
+            arizona_test_dust,
+            ip.deposition,
+            S_i,
+            T;
+            ndrange,
+        )
 
         # test if dust_activated_number_fraction is callable and returns reasonable values
         @test Array(output)[1] ≈ FT(0.0129835639)
@@ -1029,7 +1168,8 @@ function test_gpu(FT)
         dSi_dt = FT(0.03)
         N_aero = FT(3000)
 
-        kernel! = IceNucleation_MohlerDepositionRate_kernel!(backend, work_groups)
+        kernel! =
+            IceNucleation_MohlerDepositionRate_kernel!(backend, work_groups)
         kernel!(
             output,
             desert_dust,
@@ -1115,13 +1255,15 @@ function test_gpu(FT)
         T = ArrayType([FT(220)])
         Delta_a_w = ArrayType([FT(0.2907389666103033)])
 
-        kernel! = IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups)
+        kernel! =
+            IceNucleation_homogeneous_J_cubic_kernel!(backend, work_groups)
         kernel!(ip, output, Delta_a_w; ndrange)
 
         # test homogeneous_J_cubic is callable and returns a reasonable value
         @test Array(output)[1] ≈ FT(2.66194650334444e12)
 
-        kernel! = IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups)
+        kernel! =
+            IceNucleation_homogeneous_J_linear_kernel!(backend, work_groups)
         kernel!(ip, output, Delta_a_w; ndrange)
 
         # test homogeneous_J_linear is callable and returns a reasonable value
@@ -1139,7 +1281,15 @@ function test_gpu(FT)
         temp = ArrayType([FT(208), FT(208)])
 
         kernel! = h2so4_nucleation_kernel!(backend, work_groups)
-        kernel!(h2so4_nuc, output, h2so4_conc, nh3_conc, negative_ion_conc, temp; ndrange)
+        kernel!(
+            h2so4_nuc,
+            output,
+            h2so4_conc,
+            nh3_conc,
+            negative_ion_conc,
+            temp;
+            ndrange,
+        )
 
         @test all(Array(output) .> FT(0))
 
@@ -1225,10 +1375,26 @@ function test_gpu(FT)
         # test if all output is positive...
         @test all(Array(output) .> FT(0))
         #... and returns reasonable numbers
-        @test isapprox(Array(output)[1, 1], FT(0.4946323381999426 * 1e-3), rtol = 1e-2)
-        @test isapprox(Array(output)[2, 1], FT(0.26151186272014415 * 1e-3), rtol = 1e-2)
-        @test isapprox(Array(output)[1, 2], FT(1.7400778369620664 * 1e-3), rtol = 1e-2)
-        @test isapprox(Array(output)[2, 2], FT(0.11516682512848 * 1e-3), rtol = 1e-2)
+        @test isapprox(
+            Array(output)[1, 1],
+            FT(0.4946323381999426 * 1e-3),
+            rtol = 1e-2,
+        )
+        @test isapprox(
+            Array(output)[2, 1],
+            FT(0.26151186272014415 * 1e-3),
+            rtol = 1e-2,
+        )
+        @test isapprox(
+            Array(output)[1, 2],
+            FT(1.7400778369620664 * 1e-3),
+            rtol = 1e-2,
+        )
+        @test isapprox(
+            Array(output)[2, 2],
+            FT(0.11516682512848 * 1e-3),
+            rtol = 1e-2,
+        )
     end
 end
 
diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl
index c6c095674..09883369a 100644
--- a/test/heterogeneous_ice_nucleation_tests.jl
+++ b/test/heterogeneous_ice_nucleation_tests.jl
@@ -129,10 +129,12 @@ function test_heterogeneous_ice_nucleation(FT)
         for dust in [feldspar, ferrihydrite, kaolinite]
             TT.@test CMI_het.deposition_J(
                 dust,
-                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold_1) - CO.a_w_ice(tps, T_cold_1),
+                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold_1) -
+                CO.a_w_ice(tps, T_cold_1),
             ) > CMI_het.deposition_J(
                 dust,
-                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm_1) - CO.a_w_ice(tps, T_warm_1),
+                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm_1) -
+                CO.a_w_ice(tps, T_warm_1),
             )
 
             TT.@test CMI_het.deposition_J(
@@ -182,10 +184,12 @@ function test_heterogeneous_ice_nucleation(FT)
         for dust in [illite, kaolinite, desert_dust]
             TT.@test CMI_het.ABIFM_J(
                 dust,
-                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold_1) - CO.a_w_ice(tps, T_cold_1),
+                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold_1) -
+                CO.a_w_ice(tps, T_cold_1),
             ) > CMI_het.ABIFM_J(
                 dust,
-                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm_1) - CO.a_w_ice(tps, T_warm_1),
+                CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm_1) -
+                CO.a_w_ice(tps, T_warm_1),
             )
 
             TT.@test CMI_het.ABIFM_J(
@@ -225,8 +229,11 @@ function test_heterogeneous_ice_nucleation(FT)
         frequencies = FT.([0.26, 0.08])
 
         for (T, INPC, frequency) in zip(temperatures, INPCs, frequencies)
-            TT.@test CMI_het.INP_concentration_frequency(ip_frostenberg, INPC, T) ≈
-                     frequency rtol = 0.1
+            TT.@test CMI_het.INP_concentration_frequency(
+                ip_frostenberg,
+                INPC,
+                T,
+            ) ≈ frequency rtol = 0.1
         end
     end
 end
diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl
index 180c16176..2b7b6bdfc 100644
--- a/test/homogeneous_ice_nucleation_tests.jl
+++ b/test/homogeneous_ice_nucleation_tests.jl
@@ -27,10 +27,12 @@ function test_homogeneous_J_cubic(FT)
         # higher nucleation rate at colder temperatures
         TT.@test CMH.homogeneous_J_cubic(
             ip.homogeneous,
-            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) - CO.a_w_ice(tps, T_cold),
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) -
+            CO.a_w_ice(tps, T_cold),
         ) > CMH.homogeneous_J_cubic(
             ip.homogeneous,
-            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) - CO.a_w_ice(tps, T_warm),
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) -
+            CO.a_w_ice(tps, T_warm),
         )
 
         # If Δa_w out of range
@@ -60,10 +62,12 @@ function test_homogeneous_J_linear(FT)
         # higher nucleation rate at colder temperatures
         TT.@test CMH.homogeneous_J_linear(
             ip.homogeneous,
-            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) - CO.a_w_ice(tps, T_cold),
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_cold) -
+            CO.a_w_ice(tps, T_cold),
         ) > CMH.homogeneous_J_linear(
             ip.homogeneous,
-            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) - CO.a_w_ice(tps, T_warm),
+            CO.a_w_xT(H2SO4_prs, tps, x_sulph, T_warm) -
+            CO.a_w_ice(tps, T_warm),
         )
 
     end
diff --git a/test/ice_nucleation_calibration.jl b/test/ice_nucleation_calibration.jl
index ebee8dbf5..20957c79a 100644
--- a/test/ice_nucleation_calibration.jl
+++ b/test/ice_nucleation_calibration.jl
@@ -38,31 +38,49 @@ function test_J_calibration(FT, IN_mode)
 
     EKI_calibrated_parameters = EKI_output[1]
     UKI_calibrated_parameters = UKI_output[1]
-    EKI_calibrated_soln =
-        run_model([params], IN_mode, EKI_calibrated_parameters, FT, [IC], end_sim)
-    UKI_calibrated_soln =
-        run_model([params], IN_mode, UKI_calibrated_parameters, FT, [IC], end_sim)
+    EKI_calibrated_soln = run_model(
+        [params],
+        IN_mode,
+        EKI_calibrated_parameters,
+        FT,
+        [IC],
+        end_sim,
+    )
+    UKI_calibrated_soln = run_model(
+        [params],
+        IN_mode,
+        UKI_calibrated_parameters,
+        FT,
+        [IC],
+        end_sim,
+    )
     true_soln = run_model([params], IN_mode, coeff_true, FT, [IC], end_sim)
 
     TT.@testset "EKI Perfect Model Calibrations on AIDA" begin
         # test that end ICNC are similar
         if IN_mode == "ABDINM"
-            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         elseif IN_mode == "ABIFM"
-            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         elseif IN_mode == "ABHOM"
-            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test EKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         end
     end
 
     TT.@testset "UKI Perfect Model Calibrations on AIDA" begin
         # test that end ICNC are similar
         if IN_mode == "ABDINM"
-            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         elseif IN_mode == "ABIFM"
-            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         elseif IN_mode == "ABHOM"
-            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol = FT(0.3)
+            TT.@test UKI_calibrated_soln[9, end] ≈ true_soln[9, end] rtol =
+                FT(0.3)
         end
     end
 
diff --git a/test/microphysics1M_tests.jl b/test/microphysics1M_tests.jl
index 7925103a5..8e4e4eeaf 100644
--- a/test/microphysics1M_tests.jl
+++ b/test/microphysics1M_tests.jl
@@ -35,7 +35,8 @@ function test_microphysics1M(FT)
             ρ_air_ground::FT,
         ) where {FT <: Real}
             rr = q_rai / (1 - q_tot)
-            vel = FT(14.34) * ρ_air_ground^FT(0.5) * ρ^-FT(0.3654) * rr^FT(0.1346)
+            vel =
+                FT(14.34) * ρ_air_ground^FT(0.5) * ρ^-FT(0.3654) * rr^FT(0.1346)
             return vel
         end
 
@@ -62,7 +63,8 @@ function test_microphysics1M(FT)
 
         #test
         TT.@test vt_rai ≈ 3.0721397260869705 rtol = 2e-6
-        TT.@test CM1.terminal_velocity(rain, Ch2022.rain, ρ, FT(0))[1] ≈ 0 atol = eps(FT)
+        TT.@test CM1.terminal_velocity(rain, Ch2022.rain, ρ, FT(0))[1] ≈ 0 atol =
+            eps(FT)
         TT.@test v_bigger > vt_rai
     end
 
@@ -78,7 +80,8 @@ function test_microphysics1M(FT)
 
         #test
         TT.@test vt_sno ≈ 0.5151154754853068 rtol = 3e-6
-        TT.@test CM1.terminal_velocity(snow, Ch2022.large_ice, ρ, FT(0)) ≈ 0 atol = eps(FT)
+        TT.@test CM1.terminal_velocity(snow, Ch2022.large_ice, ρ, FT(0)) ≈ 0 atol =
+            eps(FT)
         TT.@test v_bigger > vt_sno
     end
 
@@ -91,8 +94,8 @@ function test_microphysics1M(FT)
         q_liq_small = FT(0.5) * q_liq_threshold
         TT.@test CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq_small) == FT(0)
 
-        TT.@test CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq_small, true) ≈ FT(0.0) atol =
-            0.15 * q_liq_threshold / τ_acnv_rai
+        TT.@test CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq_small, true) ≈
+                 FT(0.0) atol = 0.15 * q_liq_threshold / τ_acnv_rai
 
         q_liq_big = FT(1.5) * q_liq_threshold
         TT.@test CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq_big) ≈
@@ -110,21 +113,34 @@ function test_microphysics1M(FT)
         τ_acnv_sno = snow.acnv1M.τ
 
         q_ice_small = FT(0.5) * q_ice_threshold
-        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_small) == FT(0)
+        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(
+            snow.acnv1M,
+            q_ice_small,
+        ) == FT(0)
 
-        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_small, true) ≈
-                 FT(0.0) atol = FT(0.15) * q_ice_threshold / τ_acnv_sno
+        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(
+            snow.acnv1M,
+            q_ice_small,
+            true,
+        ) ≈ FT(0.0) atol = FT(0.15) * q_ice_threshold / τ_acnv_sno
 
         q_ice_big = FT(1.5) * q_ice_threshold
         TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_big) ≈
                  FT(0.5) * q_ice_threshold / τ_acnv_sno
 
-        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_big, true) ≈
-                 FT(0.5) * q_ice_threshold / τ_acnv_sno atol =
+        TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(
+            snow.acnv1M,
+            q_ice_big,
+            true,
+        ) ≈ FT(0.5) * q_ice_threshold / τ_acnv_sno atol =
             FT(0.15) * q_ice_threshold / τ_acnv_sno
 
         TT.@test CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_big) ==
-                 CM1.conv_q_ice_to_q_sno_no_supersat(snow.acnv1M, q_ice_big, false)
+                 CM1.conv_q_ice_to_q_sno_no_supersat(
+            snow.acnv1M,
+            q_ice_big,
+            false,
+        )
 
     end
 
@@ -156,13 +172,18 @@ function test_microphysics1M(FT)
         q_ice = FT(0.03) * q_vap
         q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
 
-        TT.@test CM1.conv_q_ice_to_q_sno(ice, aps, tps, q, ρ, T) ≈ FT(2.5405159487552133e-9)
+        TT.@test CM1.conv_q_ice_to_q_sno(ice, aps, tps, q, ρ, T) ≈
+                 FT(2.5405159487552133e-9)
     end
 
     TT.@testset "RainLiquidAccretion" begin
 
         # eq. 5b in [Grabowski1996](@cite)
-        function accretion_empir(q_rai::FT, q_liq::FT, q_tot::FT) where {FT <: Real}
+        function accretion_empir(
+            q_rai::FT,
+            q_liq::FT,
+            q_tot::FT,
+        ) where {FT <: Real}
             rr = q_rai / (FT(1) - q_tot)
             rl = q_liq / (FT(1) - q_tot)
             return FT(2.2) * rl * rr^FT(7 / 8)
@@ -210,17 +231,38 @@ function test_microphysics1M(FT)
         q_liq = FT(5e-4)
         q_rai = FT(5e-4)
 
-        TT.@test CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ) ≈
-                 FT(1.4150106417043544e-6)
+        TT.@test CM1.accretion(
+            liquid,
+            rain,
+            blk1mvel.rain,
+            ce,
+            q_liq,
+            q_rai,
+            ρ,
+        ) ≈ FT(1.4150106417043544e-6)
         TT.@test CM1.accretion(ice, snow, blk1mvel.snow, ce, q_ice, q_sno, ρ) ≈
                  FT(2.453070979562392e-7)
-        TT.@test CM1.accretion(liquid, snow, blk1mvel.snow, ce, q_liq, q_sno, ρ) ≈
-                 FT(2.453070979562392e-7)
+        TT.@test CM1.accretion(
+            liquid,
+            snow,
+            blk1mvel.snow,
+            ce,
+            q_liq,
+            q_sno,
+            ρ,
+        ) ≈ FT(2.453070979562392e-7)
         TT.@test CM1.accretion(ice, rain, blk1mvel.rain, ce, q_ice, q_rai, ρ) ≈
                  FT(1.768763302130443e-6)
 
-        TT.@test CM1.accretion_rain_sink(rain, ice, blk1mvel.rain, ce, q_ice, q_rai, ρ) ≈
-                 FT(3.590060148920766e-5)
+        TT.@test CM1.accretion_rain_sink(
+            rain,
+            ice,
+            blk1mvel.rain,
+            ce,
+            q_ice,
+            q_rai,
+            ρ,
+        ) ≈ FT(3.590060148920766e-5)
 
         TT.@test CM1.accretion_snow_rain(
             snow,
@@ -400,13 +442,15 @@ function test_microphysics1M(FT)
         T = FT(273.15 + 2)
         ρ = FT(1.2)
         q_sno = FT(0)
-        TT.@test CM1.snow_melt(snow, blk1mvel.snow, aps, tps, q_sno, ρ, T) ≈ FT(0)
+        TT.@test CM1.snow_melt(snow, blk1mvel.snow, aps, tps, q_sno, ρ, T) ≈
+                 FT(0)
 
         # T < T_freeze -> no snow melt
         T = FT(273.15 - 2)
         ρ = FT(1.2)
         q_sno = FT(1e-4)
-        TT.@test CM1.snow_melt(snow, blk1mvel.snow, aps, tps, q_sno, ρ, T) ≈ FT(0)
+        TT.@test CM1.snow_melt(snow, blk1mvel.snow, aps, tps, q_sno, ρ, T) ≈
+                 FT(0)
 
     end
 end
diff --git a/test/microphysics2M_tests.jl b/test/microphysics2M_tests.jl
index 3165fd4c7..3f59261f1 100644
--- a/test/microphysics2M_tests.jl
+++ b/test/microphysics2M_tests.jl
@@ -23,8 +23,13 @@ function test_microphysics2M(FT)
     VarTSc = CMP.VarTimescaleAcnv(FT)
 
     # Seifert and Beheng 2006 parameters
-    override_file =
-        joinpath(pkgdir(CM), "src", "parameters", "toml", "SB2006_limiters.toml")
+    override_file = joinpath(
+        pkgdir(CM),
+        "src",
+        "parameters",
+        "toml",
+        "SB2006_limiters.toml",
+    )
     toml_dict = CP.create_toml_dict(FT; override_file)
     SB2006 = CMP.SB2006(toml_dict)
     SB2006_no_limiters = CMP.SB2006(toml_dict, false)
@@ -80,9 +85,11 @@ function test_microphysics2M(FT)
         TT.@test CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ, N_d, true) ≈
                  CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ, N_d, false) rtol = 0.2
         TT.@test CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ, N_d, true) ≈
-                 CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ, N_d, false) rtol = 0.2
+                 CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ, N_d, false) rtol =
+            0.2
         TT.@test CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ, N_d, true) ≈
-                 CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ, N_d, false) rtol = 0.2
+                 CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ, N_d, false) rtol =
+            0.2
 
     end
 
@@ -94,13 +101,23 @@ function test_microphysics2M(FT)
 
         # compare with Wood 2005 Fig 1 panel a
         function compare(scheme, input, output; eps = 0.1)
-            TT.@test CM2.conv_q_liq_to_q_rai(scheme, input * FT(1e-3), ρ, N_d) ≈ output atol =
-                eps * output
+            TT.@test CM2.conv_q_liq_to_q_rai(scheme, input * FT(1e-3), ρ, N_d) ≈
+                     output atol = eps * output
         end
         compare(KK2000, FT(0.03138461538461537), FT(2.636846054348105e-12))
         compare(KK2000, FT(0.8738461538461537), FT(9.491665962977648e-9))
-        compare(B1994, FT(0.13999999999999999), FT(4.584323122458155e-12), eps = 1)
-        compare(B1994, FT(0.9000000000000006), FT(5.4940586176564715e-8), eps = 1)
+        compare(
+            B1994,
+            FT(0.13999999999999999),
+            FT(4.584323122458155e-12),
+            eps = 1,
+        )
+        compare(
+            B1994,
+            FT(0.9000000000000006),
+            FT(5.4940586176564715e-8),
+            eps = 1,
+        )
         compare(TC1980, FT(0.2700000000000001), FT(3.2768635256661366e-8))
         compare(TC1980, FT(0.9000000000000006), FT(5.340418612468997e-7))
         compare(LD2004, FT(0.3700000000000002), FT(8.697439193234471e-9))
@@ -108,8 +125,12 @@ function test_microphysics2M(FT)
 
         # compare with Wood 2005 Fig 1 panel b
         function compare_Nd(scheme, input, output; eps = 0.1)
-            TT.@test CM2.conv_q_liq_to_q_rai(scheme, q_liq, ρ, input * FT(1e6)) ≈ output atol =
-                eps * output
+            TT.@test CM2.conv_q_liq_to_q_rai(
+                scheme,
+                q_liq,
+                ρ,
+                input * FT(1e6),
+            ) ≈ output atol = eps * output
         end
         compare_Nd(KK2000, FT(16.13564081404141), FT(6.457285532394289e-8))
         compare_Nd(KK2000, FT(652.093931356625), FT(8.604011482409198e-11))
@@ -118,7 +139,12 @@ function test_microphysics2M(FT)
         compare_Nd(TC1980, FT(13.658073017575544), FT(2.7110779872658386e-7))
         compare_Nd(TC1980, FT(205.0970632305975), FT(1.0928660431622176e-7))
         compare_Nd(LD2004, FT(15.122629721719655), FT(1.1647783461546477e-7))
-        compare_Nd(LD2004, FT(149.01220754857331), FT(1.3917890403908125e-8), eps = 1)
+        compare_Nd(
+            LD2004,
+            FT(149.01220754857331),
+            FT(1.3917890403908125e-8),
+            eps = 1,
+        )
 
     end
 
@@ -157,9 +183,20 @@ function test_microphysics2M(FT)
 
             #action
             au = CM2.autoconversion(SB.acnv, SB.pdf_c, q_liq, q_rai, ρ, N_liq)
-            sc = CM2.liquid_self_collection(SB.acnv, SB.pdf_c, q_liq, ρ, au.dN_liq_dt)
-            au_sc =
-                CM2.autoconversion_and_liquid_self_collection(SB, q_liq, q_rai, ρ, N_liq)
+            sc = CM2.liquid_self_collection(
+                SB.acnv,
+                SB.pdf_c,
+                q_liq,
+                ρ,
+                au.dN_liq_dt,
+            )
+            au_sc = CM2.autoconversion_and_liquid_self_collection(
+                SB,
+                q_liq,
+                q_rai,
+                ρ,
+                N_liq,
+            )
 
             Lc = ρ * q_liq
             Lr = ρ * q_rai
@@ -175,7 +212,8 @@ function test_microphysics2M(FT)
             dqcdt_au = -dqrdt_au
             dNcdt_au = 2 / x_star * ρ * dqcdt_au
             dNrdt_au = -0.5 * dNcdt_au
-            dNcdt_sc = -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * Lc^2 - au.dN_liq_dt
+            dNcdt_sc =
+                -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * Lc^2 - au.dN_liq_dt
 
             #test
             TT.@test au isa CM2.LiqRaiRates
@@ -193,9 +231,20 @@ function test_microphysics2M(FT)
 
             #action
             au = CM2.autoconversion(SB.acnv, SB.pdf_c, FT(0), FT(0), ρ, N_liq)
-            sc = CM2.liquid_self_collection(SB.acnv, SB.pdf_c, FT(0), ρ, au.dN_liq_dt)
-            au_sc =
-                CM2.autoconversion_and_liquid_self_collection(SB, FT(0), FT(0), ρ, N_liq)
+            sc = CM2.liquid_self_collection(
+                SB.acnv,
+                SB.pdf_c,
+                FT(0),
+                ρ,
+                au.dN_liq_dt,
+            )
+            au_sc = CM2.autoconversion_and_liquid_self_collection(
+                SB,
+                FT(0),
+                FT(0),
+                ρ,
+                N_liq,
+            )
 
             #test
             TT.@test au.dq_liq_dt ≈ FT(0) atol = eps(FT)
@@ -266,18 +315,22 @@ function test_microphysics2M(FT)
             ρ0 = SB.pdf_r.ρ0
 
             #action
-            sc_rai = CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai)
-            br_rai = CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai)
-            sc_br_rai = CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai)
+            sc_rai =
+                CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai)
+            br_rai =
+                CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai)
+            sc_br_rai =
+                CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai)
 
             λr = CM2.pdf_rain_parameters(SB.pdf_r, q_rai, ρ, N_rai).Br
 
-            dNrdt_sc = -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ)
+            dNrdt_sc =
+                -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ)
 
             Dr =
                 (
-                    CM2.pdf_rain_parameters(SB.pdf_r, q_rai, ρ, N_rai).xr / 1000 / FT(π) *
-                    6
+                    CM2.pdf_rain_parameters(SB.pdf_r, q_rai, ρ, N_rai).xr /
+                    1000 / FT(π) * 6
                 )^FT(1 / 3)
             ΔDr = Dr - Deq
             ϕ_br =
@@ -288,8 +341,13 @@ function test_microphysics2M(FT)
 
             #test
             TT.@test sc_rai ≈ dNrdt_sc rtol = 1e-6
-            TT.@test CM2.rain_self_collection(SB.pdf_r, SB.self, FT(0), ρ, N_rai) ≈ FT(0) atol =
-                eps(FT)
+            TT.@test CM2.rain_self_collection(
+                SB.pdf_r,
+                SB.self,
+                FT(0),
+                ρ,
+                N_rai,
+            ) ≈ FT(0) atol = eps(FT)
             TT.@test br_rai ≈ dNrdt_br rtol = 1e-6
             TT.@test sc_br_rai isa NamedTuple
             TT.@test sc_br_rai.sc ≈ dNrdt_sc rtol = 1e-6
@@ -299,9 +357,12 @@ function test_microphysics2M(FT)
             q_rai = FT(0)
 
             #action
-            sc_rai = CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai)
-            br_rai = CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai)
-            sc_br_rai = CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai)
+            sc_rai =
+                CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai)
+            br_rai =
+                CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai)
+            sc_br_rai =
+                CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai)
 
             #test
             TT.@test sc_rai ≈ FT(0) atol = eps(FT)
@@ -331,10 +392,20 @@ function test_microphysics2M(FT)
         TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6
         TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6
 
-        TT.@test CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rai, ρ, FT(0))[1] ≈ 0 atol =
-            eps(FT)
-        TT.@test CM2.rain_terminal_velocity(SB2006, SB2006Vel, FT(0), ρ, N_rai)[2] ≈ 0 atol =
-            eps(FT)
+        TT.@test CM2.rain_terminal_velocity(
+            SB2006,
+            SB2006Vel,
+            q_rai,
+            ρ,
+            FT(0),
+        )[1] ≈ 0 atol = eps(FT)
+        TT.@test CM2.rain_terminal_velocity(
+            SB2006,
+            SB2006Vel,
+            FT(0),
+            ρ,
+            N_rai,
+        )[2] ≈ 0 atol = eps(FT)
     end
 
     TT.@testset "2M_microphysics - Seifert and Beheng 2006 modified rain terminal velocity without limiters" begin
@@ -346,9 +417,21 @@ function test_microphysics2M(FT)
         (; ρ0, aR, bR, cR) = SB2006Vel
 
         #action
-        vt_rai = CM2.rain_terminal_velocity(SB2006_no_limiters, SB2006Vel, q_rai, ρ, N_rai)
-
-        λr = CM2.pdf_rain_parameters(SB2006_no_limiters.pdf_r, q_rai, ρ, N_rai).λr
+        vt_rai = CM2.rain_terminal_velocity(
+            SB2006_no_limiters,
+            SB2006Vel,
+            q_rai,
+            ρ,
+            N_rai,
+        )
+
+        λr =
+            CM2.pdf_rain_parameters(
+                SB2006_no_limiters.pdf_r,
+                q_rai,
+                ρ,
+                N_rai,
+            ).λr
         _rc = -1 / (2 * cR) * log(aR / bR)
         _Γ_1(t) = exp(-t)
         _Γ_4(t) = (t^3 + 3 * t^2 + 6 * t + 6) * exp(-t)
@@ -388,18 +471,30 @@ function test_microphysics2M(FT)
 
         for SB in [SB2006, SB2006_no_limiters]
             #action
-            vt_rai = CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai, ρ, N_rai)
-            v_bigger = CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai * 2, ρ, N_rai)
+            vt_rai =
+                CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai, ρ, N_rai)
+            v_bigger =
+                CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai * 2, ρ, N_rai)
 
             #test
             TT.@test vt_rai isa Tuple
             TT.@test vt_rai[1] ≈ 1.0738503635546666 rtol = 1e-6
             TT.@test vt_rai[2] ≈ 4.00592218028957 rtol = 1e-6
 
-            TT.@test CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai, ρ, FT(0))[1] ≈ 0 atol =
-                eps(FT)
-            TT.@test CM2.rain_terminal_velocity(SB, Chen2022Vel, FT(0), ρ, N_rai)[2] ≈ 0 atol =
-                eps(FT)
+            TT.@test CM2.rain_terminal_velocity(
+                SB,
+                Chen2022Vel,
+                q_rai,
+                ρ,
+                FT(0),
+            )[1] ≈ 0 atol = eps(FT)
+            TT.@test CM2.rain_terminal_velocity(
+                SB,
+                Chen2022Vel,
+                FT(0),
+                ρ,
+                N_rai,
+            )[2] ≈ 0 atol = eps(FT)
 
             TT.@test v_bigger[1] > vt_rai[1]
             TT.@test v_bigger[2] > vt_rai[2]
@@ -444,10 +539,26 @@ function test_microphysics2M(FT)
             TT.@test evap isa NamedTuple
             TT.@test evap.evap_rate_0 ≈ evap0 rtol = 1e-4
             TT.@test evap.evap_rate_1 ≈ evap1 rtol = 1e-5
-            TT.@test CM2.rain_evaporation(SB, aps, tps, q, q_rai, ρ, FT(0), T).evap_rate_0 ≈
-                     0 atol = eps(FT)
-            TT.@test CM2.rain_evaporation(SB, aps, tps, q, FT(0), ρ, N_rai, T).evap_rate_1 ≈
-                     0 atol = eps(FT)
+            TT.@test CM2.rain_evaporation(
+                SB,
+                aps,
+                tps,
+                q,
+                q_rai,
+                ρ,
+                FT(0),
+                T,
+            ).evap_rate_0 ≈ 0 atol = eps(FT)
+            TT.@test CM2.rain_evaporation(
+                SB,
+                aps,
+                tps,
+                q,
+                FT(0),
+                ρ,
+                N_rai,
+                T,
+            ).evap_rate_1 ≈ 0 atol = eps(FT)
         end
 
         # test limit case: xr = 0 for SB with no limiters
@@ -480,7 +591,8 @@ function test_microphysics2M(FT)
         λr_l = CM2.pdf_rain_parameters(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).λr
         N₀r_l = CM2.pdf_rain_parameters(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).N₀r
         # not limited distribution parameters for rain
-        (; λr, N₀r, Ar, Br) = CM2.pdf_rain_parameters(SB2006_no_limiters.pdf_r, qᵣ, ρₐ, Nᵣ)
+        (; λr, N₀r, Ar, Br) =
+            CM2.pdf_rain_parameters(SB2006_no_limiters.pdf_r, qᵣ, ρₐ, Nᵣ)
 
         # mass of liquid droplet as a function of its diameter
         m(D) = FT(π / 6) * ρₗ * D^3
@@ -503,9 +615,15 @@ function test_microphysics2M(FT)
         m₀ = m(D₀)
         m∞ = m(D∞)
         # integral bounds computed based on the size distribution
-        D_max_limited = CM2.get_size_distribution_bound(SB2006.pdf_r, qᵣ, Nᵣ, ρₐ, eps(FT))
-        D_max =
-            CM2.get_size_distribution_bound(SB2006_no_limiters.pdf_r, qᵣ, Nᵣ, ρₐ, eps(FT))
+        D_max_limited =
+            CM2.get_size_distribution_bound(SB2006.pdf_r, qᵣ, Nᵣ, ρₐ, eps(FT))
+        D_max = CM2.get_size_distribution_bound(
+            SB2006_no_limiters.pdf_r,
+            qᵣ,
+            Nᵣ,
+            ρₐ,
+            eps(FT),
+        )
 
         # Sanity checks for number concentrations for rain
         ND = QGK.quadgk(x -> f_D(x), D₀, D∞)[1]
@@ -513,7 +631,13 @@ function test_microphysics2M(FT)
         ND_lim = QGK.quadgk(x -> f_D_limited(x), D₀, D∞)[1]
         Nx_lim = QGK.quadgk(x -> f_x_limited(x), m₀, m∞)[1]
         ND_bounded = QGK.quadgk(
-            x -> CM2.size_distribution(SB2006_no_limiters.pdf_r, FT(x), qᵣ, ρₐ, Nᵣ),
+            x -> CM2.size_distribution(
+                SB2006_no_limiters.pdf_r,
+                FT(x),
+                qᵣ,
+                ρₐ,
+                Nᵣ,
+            ),
             0,
             D_max,
         )[1]
@@ -542,14 +666,21 @@ function test_microphysics2M(FT)
         qD_bounded =
             QGK.quadgk(
                 x ->
-                    m(x) *
-                    CM2.size_distribution(SB2006_no_limiters.pdf_r, FT(x), qᵣ, ρₐ, Nᵣ),
+                    m(x) * CM2.size_distribution(
+                        SB2006_no_limiters.pdf_r,
+                        FT(x),
+                        qᵣ,
+                        ρₐ,
+                        Nᵣ,
+                    ),
                 0,
                 D_max,
             )[1] / ρₐ
         qD_bounded_limited =
             QGK.quadgk(
-                x -> m(x) * CM2.size_distribution(SB2006.pdf_r, FT(x), qᵣ, ρₐ, Nᵣ),
+                x ->
+                    m(x) *
+                    CM2.size_distribution(SB2006.pdf_r, FT(x), qᵣ, ρₐ, Nᵣ),
                 0,
                 D_max_limited,
             )[1] / ρₐ
@@ -618,8 +749,13 @@ function test_microphysics2M(FT)
         TT.@test qD ≈ qₗ rtol = FT(1e-6)
 
         # integral bounds computed based on the size distribution
-        D_max =
-            CM2.get_size_distribution_bound(SB2006_no_limiters.pdf_c, qₗ, Nₗ, ρₐ, eps(FT))
+        D_max = CM2.get_size_distribution_bound(
+            SB2006_no_limiters.pdf_c,
+            qₗ,
+            Nₗ,
+            ρₐ,
+            eps(FT),
+        )
         # Sanity checks specific humidity and number concentration with diameter distribution
         qD_bounded =
             QGK.quadgk(
@@ -627,12 +763,24 @@ function test_microphysics2M(FT)
                     FT(π / 6) *
                     ρₗ *
                     x^3 *
-                    CM2.size_distribution(SB2006_no_limiters.pdf_c, FT(x), qₗ, ρₐ, Nₗ),
+                    CM2.size_distribution(
+                        SB2006_no_limiters.pdf_c,
+                        FT(x),
+                        qₗ,
+                        ρₐ,
+                        Nₗ,
+                    ),
                 0,
                 D_max,
             )[1] / ρₐ
         ND_bounded = QGK.quadgk(
-            x -> CM2.size_distribution(SB2006_no_limiters.pdf_c, FT(x), qₗ, ρₐ, Nₗ),
+            x -> CM2.size_distribution(
+                SB2006_no_limiters.pdf_c,
+                FT(x),
+                qₗ,
+                ρₐ,
+                Nₗ,
+            ),
             0,
             D_max,
         )[1]
diff --git a/test/p3_tests.jl b/test/p3_tests.jl
index 354a2b192..6b1827f4e 100644
--- a/test/p3_tests.jl
+++ b/test/p3_tests.jl
@@ -25,7 +25,8 @@ function test_incomplete_gamma_function_approximation(FT)
 
         for it in [1, 2, 3]
             for jt in [1, 2, 3]
-                TT.@test P3.Γ_lower(a[it], z[jt]) ≈ Γ_lower_reference[it][jt] rtol = 1e-3
+                TT.@test P3.Γ_lower(a[it], z[jt]) ≈ Γ_lower_reference[it][jt] rtol =
+                    1e-3
             end
         end
     end
@@ -45,7 +46,11 @@ function test_thresholds_solver(FT)
 
         # test asserts
         for _ρ_r in (FT(0), FT(-1))
-            TT.@test_throws AssertionError("ρ_r > FT(0)") P3.thresholds(p3, _ρ_r, F_rim)
+            TT.@test_throws AssertionError("ρ_r > FT(0)") P3.thresholds(
+                p3,
+                _ρ_r,
+                F_rim,
+            )
         end
         for _ρ_r in (FT(0), FT(-1))
             TT.@test_throws AssertionError("ρ_r <= p3.ρ_l") P3.thresholds(
@@ -56,10 +61,18 @@ function test_thresholds_solver(FT)
         end
 
         for _F_rim in (FT(-1 * eps(FT)), FT(-1))
-            TT.@test_throws AssertionError("F_rim >= FT(0)") P3.thresholds(p3, ρ_r, _F_rim)
+            TT.@test_throws AssertionError("F_rim >= FT(0)") P3.thresholds(
+                p3,
+                ρ_r,
+                _F_rim,
+            )
         end
         for _F_rim in (FT(1), FT(1.5))
-            TT.@test_throws AssertionError("F_rim < FT(1)") P3.thresholds(p3, ρ_r, _F_rim)
+            TT.@test_throws AssertionError("F_rim < FT(1)") P3.thresholds(
+                p3,
+                ρ_r,
+                _F_rim,
+            )
         end
 
         # Test if the P3 scheme solution satisifies the conditions
@@ -68,7 +81,8 @@ function test_thresholds_solver(FT)
             for ρ_r in ρ_r_good
                 sol = P3.thresholds(p3, ρ_r, F_rim)
                 atol = 5e3 * eps(FT)
-                TT.@test sol.D_cr ≈ P3.D_cr_helper(p3, F_rim, sol[3]) atol = atol
+                TT.@test sol.D_cr ≈ P3.D_cr_helper(p3, F_rim, sol[3]) atol =
+                    atol
                 TT.@test sol.D_gr ≈ P3.D_gr_helper(p3, sol[3]) atol = atol
                 TT.@test sol.ρ_g ≈ P3.ρ_g_helper(ρ_r, F_rim, sol[4]) atol = atol
                 TT.@test sol.ρ_d ≈ P3.ρ_d_helper(p3, sol[1], sol[2]) atol = atol
@@ -117,13 +131,15 @@ function test_thresholds_solver(FT)
         TT.@test P3.p3_area(p3, D_1, F_rim, F_liq, th) == P3.A_s(D_1)
         TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) == P3.A_ns(p3, D_2)
         TT.@test P3.p3_area(p3, D_3, F_rim, F_liq, th) == P3.A_s(D_3)
-        TT.@test P3.p3_area(p3, D_cr, F_rim, F_liq, th) == P3.A_r(p3, F_rim, D_cr)
+        TT.@test P3.p3_area(p3, D_cr, F_rim, F_liq, th) ==
+                 P3.A_r(p3, F_rim, D_cr)
 
         # test mass
         TT.@test P3.p3_mass(p3, D_1, F_rim, F_liq, th) == P3.mass_s(D_1, p3.ρ_i)
         TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) == P3.mass_nl(p3, D_2)
         TT.@test P3.p3_mass(p3, D_3, F_rim, F_liq, th) == P3.mass_s(D_3, th.ρ_g)
-        TT.@test P3.p3_mass(p3, D_cr, F_rim, F_liq, th) == P3.mass_r(p3, D_cr, F_rim)
+        TT.@test P3.p3_mass(p3, D_cr, F_rim, F_liq, th) ==
+                 P3.mass_r(p3, D_cr, F_rim)
 
         # test density
         TT.@test P3.p3_density(p3, D_1, F_rim, th) ≈ p3.ρ_i
@@ -170,20 +186,25 @@ function test_thresholds_solver(FT)
 
         # test mass
         TT.@test P3.p3_mass(p3, D_1, F_rim, F_liq, th) ==
-                 (1 - F_liq) * P3.mass_s(D_1, p3.ρ_i) + F_liq * P3.mass_s(D_1, p3.ρ_l)
+                 (1 - F_liq) * P3.mass_s(D_1, p3.ρ_i) +
+                 F_liq * P3.mass_s(D_1, p3.ρ_l)
         TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) ==
-                 (1 - F_liq) * P3.mass_nl(p3, D_2) + F_liq * P3.mass_s(D_2, p3.ρ_l)
+                 (1 - F_liq) * P3.mass_nl(p3, D_2) +
+                 F_liq * P3.mass_s(D_2, p3.ρ_l)
         TT.@test P3.p3_mass(p3, D_3, F_rim, F_liq, th) ==
-                 (1 - F_liq) * P3.mass_s(D_3, th.ρ_g) + F_liq * P3.mass_s(D_3, p3.ρ_l)
+                 (1 - F_liq) * P3.mass_s(D_3, th.ρ_g) +
+                 F_liq * P3.mass_s(D_3, p3.ρ_l)
         TT.@test P3.p3_mass(p3, D_cr, F_rim, F_liq, th) ==
-                 (1 - F_liq) * P3.mass_r(p3, D_cr, F_rim) + F_liq * P3.mass_s(D_cr, p3.ρ_l)
+                 (1 - F_liq) * P3.mass_r(p3, D_cr, F_rim) +
+                 F_liq * P3.mass_s(D_cr, p3.ρ_l)
 
         # test F_rim = 0 and D > D_th
         F_rim = FT(0)
         TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) ==
                  (1 - F_liq) * P3.A_ns(p3, D_2) + F_liq * P3.A_s(D_2)
         TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) ==
-                 (1 - F_liq) * P3.mass_nl(p3, D_2) + F_liq * P3.mass_s(D_2, p3.ρ_l)
+                 (1 - F_liq) * P3.mass_nl(p3, D_2) +
+                 F_liq * P3.mass_s(D_2, p3.ρ_l)
     end
 end
 
@@ -216,7 +237,9 @@ function test_shape_solver(FT)
                             # Convert λ to ensure it remains positive
                             x = log(λ_ex)
                             # Compute mass density based on input shape parameters
-                            L_calc = N * P3.L_over_N_gamma(p3, F_liq, F_rim, x, μ_ex, th)
+                            L_calc =
+                                N *
+                                P3.L_over_N_gamma(p3, F_liq, F_rim, x, μ_ex, th)
 
                             if L_calc < FT(1)
                                 # Solve for shape parameters
@@ -371,10 +394,28 @@ function test_bulk_terminal_velocities(FT)
         for k in 1:length(F_rims)
             F_rim = F_rims[k]
             F_liq = FT(0)
-            vel =
-                P3.ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_rim, F_liq, ρ_a, false)
-            vel_ϕ =
-                P3.ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_rim, F_liq, ρ_a, true)
+            vel = P3.ice_terminal_velocity(
+                p3,
+                Chen2022,
+                L,
+                N,
+                ρ_r,
+                F_rim,
+                F_liq,
+                ρ_a,
+                false,
+            )
+            vel_ϕ = P3.ice_terminal_velocity(
+                p3,
+                Chen2022,
+                L,
+                N,
+                ρ_r,
+                F_rim,
+                F_liq,
+                ρ_a,
+                true,
+            )
 
             # number weighted
             TT.@test vel[1] > 0
@@ -404,10 +445,28 @@ function test_bulk_terminal_velocities(FT)
         for k in 1:length(F_liqs)
             F_liq = F_liqs[k]
             F_rim = FT(0.4)
-            vel =
-                P3.ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_rim, F_liq, ρ_a, false)
-            vel_ϕ =
-                P3.ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_rim, F_liq, ρ_a, true)
+            vel = P3.ice_terminal_velocity(
+                p3,
+                Chen2022,
+                L,
+                N,
+                ρ_r,
+                F_rim,
+                F_liq,
+                ρ_a,
+                false,
+            )
+            vel_ϕ = P3.ice_terminal_velocity(
+                p3,
+                Chen2022,
+                L,
+                N,
+                ρ_r,
+                F_rim,
+                F_liq,
+                ρ_a,
+                true,
+            )
             # number weighted
             TT.@test vel[1] > 0
             TT.@test vel_ϕ[1] > 0
@@ -464,7 +523,14 @@ function test_numerical_integrals(FT)
                 L = Ls[i]
 
                 # Get shape parameters, thresholds and intergal bounds
-                λ, N_0 = P3.distribution_parameter_solver(p3, L, N, ρ_r, F_rim, F_liq)
+                λ, N_0 = P3.distribution_parameter_solver(
+                    p3,
+                    L,
+                    N,
+                    ρ_r,
+                    F_rim,
+                    F_liq,
+                )
                 th = P3.thresholds(p3, ρ_r, F_rim)
                 ice_bound = P3.get_ice_bound(p3, λ, N, tolerance)
 
@@ -508,7 +574,10 @@ function test_numerical_integrals(FT)
 
                 # Dₘ comparisons
                 D_m = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq)
-                h(d) = d * P3.p3_mass(p3, d, F_rim, F_liq, th) * P3.N′ice(p3, d, λ, N_0)
+                h(d) =
+                    d *
+                    P3.p3_mass(p3, d, F_rim, F_liq, th) *
+                    P3.N′ice(p3, d, λ, N_0)
                 qgk_D_m, = QGK.quadgk(d -> h(d) / L, FT(0), ice_bound)
                 TT.@test D_m ≈ qgk_D_m rtol = 1e-2
             end
@@ -706,13 +775,15 @@ function test_p3_melting(FT)
         F_liq = FT(0) # process not dependent on F_liq
 
         T_cold = FT(273.15 - 0.01)
-        rate = P3.ice_melt(p3, vel, aps, tps, Lᵢ, Nᵢ, T_cold, ρₐ, F_rim, ρ_rim, dt)
+        rate =
+            P3.ice_melt(p3, vel, aps, tps, Lᵢ, Nᵢ, T_cold, ρₐ, F_rim, ρ_rim, dt)
 
         TT.@test rate.dNdt == 0
         TT.@test rate.dLdt == 0
 
         T_warm = FT(273.15 + 0.01)
-        rate = P3.ice_melt(p3, vel, aps, tps, Lᵢ, Nᵢ, T_warm, ρₐ, F_rim, ρ_rim, dt)
+        rate =
+            P3.ice_melt(p3, vel, aps, tps, Lᵢ, Nᵢ, T_warm, ρₐ, F_rim, ρ_rim, dt)
 
         TT.@test rate.dNdt >= 0
         TT.@test rate.dLdt >= 0
@@ -721,7 +792,19 @@ function test_p3_melting(FT)
         TT.@test rate.dLdt ≈ 8.601633922081113e-5 rtol = 1e-6
 
         T_vwarm = FT(273.15 + 0.1)
-        rate = P3.ice_melt(p3, vel, aps, tps, Lᵢ, Nᵢ, T_vwarm, ρₐ, F_rim, ρ_rim, dt)
+        rate = P3.ice_melt(
+            p3,
+            vel,
+            aps,
+            tps,
+            Lᵢ,
+            Nᵢ,
+            T_vwarm,
+            ρₐ,
+            F_rim,
+            ρ_rim,
+            dt,
+        )
 
         TT.@test rate.dNdt == Nᵢ
         TT.@test rate.dLdt == Lᵢ
diff --git a/test/performance_tests.jl b/test/performance_tests.jl
index 9137c3e20..a63ed0307 100644
--- a/test/performance_tests.jl
+++ b/test/performance_tests.jl
@@ -23,7 +23,13 @@ import CloudMicrophysics.CloudDiagnostics as CMP
 
 @info "Performance Tests"
 
-function bench_press(foo, args, min_run_time, min_memory = 0.0, min_allocs = 0.0)
+function bench_press(
+    foo,
+    args,
+    min_run_time,
+    min_memory = 0.0,
+    min_allocs = 0.0,
+)
     println("Testing ", "$foo")
 
     # converts a string argument to tuple
@@ -47,7 +53,13 @@ end
 function benchmark_test(FT)
 
     # Artifact calling
-    bench_press(AFC.AIDA_ice_nucleation, ("in05_17_aida.edf"), 50000, 30000, 300)
+    bench_press(
+        AFC.AIDA_ice_nucleation,
+        ("in05_17_aida.edf"),
+        50000,
+        30000,
+        300,
+    )
 
     # 0-moment microphysics
     p0m = CMP.Parameters0M(FT)
@@ -58,8 +70,13 @@ function benchmark_test(FT)
     snow = CMP.Snow(FT)
     ce = CMP.CollisionEff(FT)
     # 2-moment microphysics
-    override_file =
-        joinpath(pkgdir(CM), "src", "parameters", "toml", "SB2006_limiters.toml")
+    override_file = joinpath(
+        pkgdir(CM),
+        "src",
+        "parameters",
+        "toml",
+        "SB2006_limiters.toml",
+    )
     toml_dict = CP.create_toml_dict(FT; override_file)
     sb2006 = CMP.SB2006(toml_dict)
     sb2006_no_limiters = CMP.SB2006(toml_dict, false)
@@ -118,8 +135,15 @@ function benchmark_test(FT)
     N_aer = FT(100.0 * 1e6)
     M_seasalt = FT(0.058443)
     κ_seasalt = FT(1.12)
-    seasalt_mode =
-        AM.Mode_κ(r_aer, σ_aer, N_aer, (FT(1),), (FT(1),), (M_seasalt,), (κ_seasalt,))
+    seasalt_mode = AM.Mode_κ(
+        r_aer,
+        σ_aer,
+        N_aer,
+        (FT(1),),
+        (FT(1),),
+        (M_seasalt,),
+        (κ_seasalt,),
+    )
     aer_distr = AM.AerosolDistribution((seasalt_mode,))
 
     x_sulph = FT(0.1)
@@ -183,9 +207,17 @@ function benchmark_test(FT)
 
     # Chen 2022 terminal velocity
     bench_press(CMN.terminal_velocity, (liquid, ch2022.rain, ρ_air, q_liq), 350)
-    bench_press(CMN.terminal_velocity, (ice, ch2022.small_ice, ρ_air, q_ice), 400)
+    bench_press(
+        CMN.terminal_velocity,
+        (ice, ch2022.small_ice, ρ_air, q_ice),
+        400,
+    )
     bench_press(CM1.terminal_velocity, (rain, ch2022.rain, ρ_air, q_rai), 750)
-    bench_press(CM1.terminal_velocity, (snow, ch2022.large_ice, ρ_air, q_sno), 750)
+    bench_press(
+        CM1.terminal_velocity,
+        (snow, ch2022.large_ice, ρ_air, q_sno),
+        750,
+    )
 
     # aerosol activation
     bench_press(
@@ -246,7 +278,11 @@ function benchmark_test(FT)
     bench_press(CM0.remove_precipitation, (p0m, q), 12)
 
     # 1-moment
-    bench_press(CM1.accretion, (liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air), 350)
+    bench_press(
+        CM1.accretion,
+        (liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air),
+        350,
+    )
     bench_press(CMD.radar_reflectivity_1M, (rain, q_rai, ρ_air), 300)
 
     # 2-moment
@@ -256,13 +292,21 @@ function benchmark_test(FT)
             (sb, q_liq, q_rai, ρ_air, N_liq),
             300,
         )
-        bench_press(CM2.rain_self_collection_and_breakup, (sb, q_rai, ρ_air, N_rai), 1200)
+        bench_press(
+            CM2.rain_self_collection_and_breakup,
+            (sb, q_rai, ρ_air, N_rai),
+            1200,
+        )
         bench_press(
             CM2.rain_evaporation,
             (sb, aps, tps, q, q_rai, ρ_air, N_rai, T_air),
             2000,
         )
-        bench_press(CM2.rain_terminal_velocity, (sb, sb2006vel, q_rai, ρ_air, N_rai), 700)
+        bench_press(
+            CM2.rain_terminal_velocity,
+            (sb, sb2006vel, q_rai, ρ_air, N_rai),
+            700,
+        )
         bench_press(
             CM2.rain_terminal_velocity,
             (sb, ch2022.rain, q_rai, ρ_air, N_rai),
@@ -273,7 +317,11 @@ function benchmark_test(FT)
             (sb, q_liq, q_rai, N_liq, N_rai, ρ_air),
             2000,
         )
-        bench_press(CMD.effective_radius_2M, (sb, q_liq, q_rai, N_liq, N_rai, ρ_air), 2000)
+        bench_press(
+            CMD.effective_radius_2M,
+            (sb, q_liq, q_rai, N_liq, N_rai, ρ_air),
+            2000,
+        )
     end
     bench_press(
         CMD.effective_radius_Liu_Hallet_97,
@@ -282,7 +330,11 @@ function benchmark_test(FT)
     )
     # Homogeneous Nucleation
     bench_press(HN.h2so4_nucleation_rate, (1e12, 1.0, 1.0, 208, h2so4_nuc), 470)
-    bench_press(HN.organic_nucleation_rate, (0.0, 1e3, 1e3, 1e3, 300, 1, organ_nuc), 850)
+    bench_press(
+        HN.organic_nucleation_rate,
+        (0.0, 1e3, 1e3, 1e3, 300, 1, organ_nuc),
+        850,
+    )
     bench_press(
         HN.organic_and_h2so4_nucleation_rate,
         (2.6e6, 1.0, 1.0, 300, 1, mixed_nuc),
diff --git a/test/precipitation_susceptibility_tests.jl b/test/precipitation_susceptibility_tests.jl
index deaf3f4b0..aaa12fcd8 100644
--- a/test/precipitation_susceptibility_tests.jl
+++ b/test/precipitation_susceptibility_tests.jl
@@ -10,7 +10,10 @@ Logarithmic derivative of universal function for autoconversion as described
     in Glassmeier & Lohmann, which is (1 + Phi(τ)/(1 - τ^2)) for Phi(τ)
     from SB2006.
 """
-d_ln_phi_au_d_ln_τ((; A, a, b)::CMP.AcnvSB2006{FT}, τ::FT) where {FT <: AbstractFloat} =
+d_ln_phi_au_d_ln_τ(
+    (; A, a, b)::CMP.AcnvSB2006{FT},
+    τ::FT,
+) where {FT <: AbstractFloat} =
     -(
         A *
         τ^a *
@@ -22,8 +25,10 @@ d_ln_phi_au_d_ln_τ((; A, a, b)::CMP.AcnvSB2006{FT}, τ::FT) where {FT <: Abstra
 Logarithmic derivative of universal function for accretion as described
     in Glassmeier & Lohmann
 """
-d_ln_phi_acc_d_ln_τ((; τ0, c)::CMP.AccrSB2006{FT}, τ::FT) where {FT <: AbstractFloat} =
-    (c * τ0) / (τ + τ0)
+d_ln_phi_acc_d_ln_τ(
+    (; τ0, c)::CMP.AccrSB2006{FT},
+    τ::FT,
+) where {FT <: AbstractFloat} = (c * τ0) / (τ + τ0)
 
 
 TT.@testset "precipitation_susceptibility_SB2006" begin
@@ -36,14 +41,27 @@ TT.@testset "precipitation_susceptibility_SB2006" begin
 
     τ = FT(1) - q_liq / (q_liq + q_rai)
 
-    aut_rates =
-        CMPS.precipitation_susceptibility_autoconversion(sb2006, q_liq, q_rai, ρ, N_liq)
+    aut_rates = CMPS.precipitation_susceptibility_autoconversion(
+        sb2006,
+        q_liq,
+        q_rai,
+        ρ,
+        N_liq,
+    )
 
-    acc_rates = CMPS.precipitation_susceptibility_accretion(sb2006, q_liq, q_rai, ρ, N_liq)
+    acc_rates = CMPS.precipitation_susceptibility_accretion(
+        sb2006,
+        q_liq,
+        q_rai,
+        ρ,
+        N_liq,
+    )
 
     TT.@test aut_rates.d_ln_pp_d_ln_N_liq ≈ -2
-    TT.@test aut_rates.d_ln_pp_d_ln_q_liq ≈ 4 - (1 - τ) * d_ln_phi_au_d_ln_τ(sb2006.acnv, τ)
-    TT.@test aut_rates.d_ln_pp_d_ln_q_rai ≈ (1 - τ) * d_ln_phi_au_d_ln_τ(sb2006.acnv, τ)
+    TT.@test aut_rates.d_ln_pp_d_ln_q_liq ≈
+             4 - (1 - τ) * d_ln_phi_au_d_ln_τ(sb2006.acnv, τ)
+    TT.@test aut_rates.d_ln_pp_d_ln_q_rai ≈
+             (1 - τ) * d_ln_phi_au_d_ln_τ(sb2006.acnv, τ)
 
     TT.@test acc_rates.d_ln_pp_d_ln_q_liq ≈
              1 - (1 - τ) * d_ln_phi_acc_d_ln_τ(sb2006.accr, τ)