diff --git a/box/Alpert_Knopf_2016_backward.jl b/box/Alpert_Knopf_2016_backward.jl
index 2a9e73459..9c0b1c9a6 100644
--- a/box/Alpert_Knopf_2016_backward.jl
+++ b/box/Alpert_Knopf_2016_backward.jl
@@ -95,7 +95,7 @@ J_immer = @. CMI_het.ABIFM_J(aerosol, Δa) # m^-2 s^-1
 
 #! format: off
 # Plot results
-fig = MK.Figure(resolution = (1200, 400))
+fig = MK.Figure(size = (1200, 400))
 ax1 = MK.Axis(fig[1, 1], ylabel = "J_immer [cm^-2 s^-1]", xlabel = "Temperature [K]", yscale = log10)
 ax2 = MK.Axis(fig[1, 2], ylabel = "frozen fraction",      xlabel = "Temperature [K]")
 ax3 = MK.Axis(fig[1, 3], ylabel = "total A [cm^2]",       xlabel = "Temperature [K]", yscale = log10)
diff --git a/box/Alpert_Knopf_2016_forward.jl b/box/Alpert_Knopf_2016_forward.jl
index b489e6a2e..4a19e98c4 100644
--- a/box/Alpert_Knopf_2016_forward.jl
+++ b/box/Alpert_Knopf_2016_forward.jl
@@ -66,7 +66,7 @@ ff_var = @. frozen_frac(sol_var[4, :], sol_var[3, :])
 
 #! format: off
 # Plot results
-fig = MK.Figure(resolution = (1200, 400))
+fig = MK.Figure(size = (1200, 400))
 ax1 = MK.Axis(fig[1, 1], ylabel = "J_immer [cm^-2 s^-1]", xlabel = "Temperature [K]", yscale = log10)
 ax2 = MK.Axis(fig[1, 2], ylabel = "frozen fraction",      xlabel = "Temperature [K]")
 ax3 = MK.Axis(fig[1, 3], ylabel = "total A [cm2]",        xlabel = "Temperature [K]", yscale = log10)
diff --git a/docs/src/plots/ARGplots.jl b/docs/src/plots/ARGplots.jl
index 06d2b784c..e46089001 100644
--- a/docs/src/plots/ARGplots.jl
+++ b/docs/src/plots/ARGplots.jl
@@ -1,19 +1,15 @@
-import Plots
 
-import CloudMicrophysics
-import ClimaParams
-import Thermodynamics
+import Plots as PL
 
-const PL = Plots
-const AM = CloudMicrophysics.AerosolModel
-const AA = CloudMicrophysics.AerosolActivation
-const CP = ClimaParams
-const CMP = CloudMicrophysics.Parameters
-const TD = Thermodynamics
+import Thermodynamics as TD
+
+import CloudMicrophysics.AerosolModel as AM
+import CloudMicrophysics.AerosolActivation as AA
+import CloudMicrophysics.Parameters as CMP
 
-FT = Float64
 
-tps = Thermodynamics.Parameters.ThermodynamicsParameters(FT)
+FT = Float64
+tps = TD.Parameters.ThermodynamicsParameters(FT)
 aip = CMP.AirProperties(FT)
 ap = CMP.AerosolActivationParameters(FT)
 
diff --git a/docs/src/plots/ARGplots_fig1.jl b/docs/src/plots/ARGplots_fig1.jl
index 51e709c74..4e3ce7955 100644
--- a/docs/src/plots/ARGplots_fig1.jl
+++ b/docs/src/plots/ARGplots_fig1.jl
@@ -1,17 +1,15 @@
-import Plots
 
-import CloudMicrophysics
-import ClimaParams
-import Thermodynamics
+import Plots as PL
+
+import Thermodynamics as TD
+
+import CloudMicrophysics.AerosolModel as AM
+import CloudMicrophysics.AerosolActivation as AA
+import CloudMicrophysics.Parameters as CMP
 
-const PL = Plots
-const AM = CloudMicrophysics.AerosolModel
-const AA = CloudMicrophysics.AerosolActivation
-const CMP = CloudMicrophysics.Parameters
-const TD = Thermodynamics
 
 FT = Float64
-tps = Thermodynamics.Parameters.ThermodynamicsParameters(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
 aip = CMP.AirProperties(FT)
 ap = CMP.AerosolActivationParameters(FT)
 
diff --git a/docs/src/plots/Baumgartner2022_fig5.jl b/docs/src/plots/Baumgartner2022_fig5.jl
index 02a59f391..d0876da20 100644
--- a/docs/src/plots/Baumgartner2022_fig5.jl
+++ b/docs/src/plots/Baumgartner2022_fig5.jl
@@ -1,12 +1,10 @@
 import RootSolvers as RS
 import CairoMakie as MK
 
-import CloudMicrophysics as CM
 import Thermodynamics as TD
-import ClimaParams as CP
 
-const CMO = CM.Common
-const CMP = CM.Parameters
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
 tps = TD.Parameters.ThermodynamicsParameters(FT)
@@ -54,12 +52,12 @@ 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
+    (CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T)) for
     T in T_range
 ]
 
 #! format: off
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 ax1 = MK.Axis(fig[1, 1], title = "Water Activity vs Temperature", ylabel = "a_w [-]", xlabel = "T [K]")
 MK.lines!(ax1, BG2022_fig5_T, BG2022_fig5_aw, label = "Baumgartner2022", linestyle = :dash, color = :blue)
 MK.lines!(ax1, T_range, a_w, label = "CloudMicrophysics", color = :blue)
@@ -70,7 +68,7 @@ MK.lines!(ax1, T_range, a_w_ice, label = "a_w_ice", color = :red)
 MK.axislegend(position = :lt)
 MK.save("Baumgartner2022_fig5.svg", fig)
 
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 ax1 = MK.Axis(fig[1, 1], title = "Saturated Vapor Pressure vs Temperature", ylabel = "Pressure [Pa]", xlabel = "T [K]")
 MK.lines!(ax1, T_range, [Baumgartner_p_ice(T) for T in T_range], label = "Baumgartner's p_ice", color = :blue, linestyle = :dash)
 MK.lines!(ax1, T_range, p_sat_ice, label = "CM's p_ice", color = :blue)
diff --git a/docs/src/plots/CLOUD_Nucleation_Plots.jl b/docs/src/plots/CLOUD_Nucleation_Plots.jl
index b1dc21ac2..19ef0f2f4 100644
--- a/docs/src/plots/CLOUD_Nucleation_Plots.jl
+++ b/docs/src/plots/CLOUD_Nucleation_Plots.jl
@@ -1,8 +1,7 @@
-using Plots
+import Plots
 
-import ClimaParams as CP
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Nucleation as Nucleation
+import CloudMicrophysics.Nucleation
 
 FT = Float64
 params = CMP.H2S04NucleationParameters(FT)
diff --git a/docs/src/plots/Cholette2019_fig1.jl b/docs/src/plots/Cholette2019_fig1.jl
index 16cce4136..d9db92765 100644
--- a/docs/src/plots/Cholette2019_fig1.jl
+++ b/docs/src/plots/Cholette2019_fig1.jl
@@ -1,11 +1,8 @@
-import ClimaParams
-import CloudMicrophysics as CM
+import CairoMakie as MK
+
 import CloudMicrophysics.P3Scheme as P3
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Common as CO
-import CairoMakie as Plt
-
-const PSP3 = CMP.ParametersP3
+import CloudMicrophysics.Parameters.ParametersP3 as PSP3
 
 FT = Float64
 
@@ -72,10 +69,10 @@ function fig1()
     colors = [:black, :blue, :red, :green]
     res = 50
 
-    fig = Plt.Figure(size = (1200, 400))
+    fig = MK.Figure(size = (1200, 400))
 
     #F_r = 0
-    ax1 = Plt.Axis(
+    ax1 = MK.Axis(
         fig[1:7, 1:9],
         title = "Fᵣ = 0",
         xlabel = "Dₚ (m)",
@@ -92,8 +89,8 @@ function fig1()
             P3.thresholds(p3, ρ_r, F_rs[1]),
             res,
         )
-        Plt.lines!(ax1, D_ps, V_0_ϕ, color = colors[i], label = labels[i])
-        Plt.lines!(
+        MK.lines!(ax1, D_ps, V_0_ϕ, color = colors[i], label = labels[i])
+        MK.lines!(
             ax1,
             D_ps,
             V_0,
@@ -104,7 +101,7 @@ function fig1()
     end
 
     # F_r = 1
-    ax2 = Plt.Axis(
+    ax2 = MK.Axis(
         fig[1:7, 10:18],
         title = "Fᵣ = 1",
         xlabel = "Dₚ (m)",
@@ -120,8 +117,8 @@ function fig1()
             P3.thresholds(p3, ρ_r, F_rs[2]),
             res,
         )
-        Plt.lines!(ax2, D_ps, V_1_ϕ, color = colors[i], label = labels[i])
-        Plt.lines!(
+        MK.lines!(ax2, D_ps, V_1_ϕ, color = colors[i], label = labels[i])
+        MK.lines!(
             ax2,
             D_ps,
             V_1,
@@ -131,20 +128,20 @@ function fig1()
         )
     end
 
-    Plt.Legend(
+    MK.Legend(
         fig[4, 19:22],
         [
-            Plt.LineElement(color = :black),
-            Plt.LineElement(color = :blue),
-            Plt.LineElement(color = :red),
-            Plt.LineElement(color = :green),
-            Plt.LineElement(color = :black, linestyle = :dash),
+            MK.LineElement(color = :black),
+            MK.LineElement(color = :blue),
+            MK.LineElement(color = :red),
+            MK.LineElement(color = :green),
+            MK.LineElement(color = :black, linestyle = :dash),
         ],
         labels,
         framevisible = false,
     )
 
-    Plt.linkaxes!(ax1, ax2)
+    MK.linkaxes!(ax1, ax2)
     ax1.xticks = (
         [0, 0.002, 0.004, 0.006, 0.008, 0.01],
         string.([0, 0.002, 0.004, 0.006, 0.008, 0.01]),
@@ -155,8 +152,8 @@ function fig1()
     )
 
 
-    Plt.resize_to_layout!(fig)
-    Plt.save("Choletteetal2019_fig1.svg", fig)
+    MK.resize_to_layout!(fig)
+    MK.save("Choletteetal2019_fig1.svg", fig)
 end
 
 fig1()
diff --git a/docs/src/plots/CloudDiagnostics.jl b/docs/src/plots/CloudDiagnostics.jl
index f952335f0..891291521 100644
--- a/docs/src/plots/CloudDiagnostics.jl
+++ b/docs/src/plots/CloudDiagnostics.jl
@@ -1,5 +1,5 @@
-import CairoMakie as PL
-PL.activate!(type = "svg")
+import CairoMakie as MK
+MK.activate!(type = "svg")
 
 import ClimaParams as CP
 
@@ -16,11 +16,7 @@ cloud_liquid = CMP.CloudLiquid(FT)
 cloud_ice = CMP.CloudIce(FT)
 
 override_file = joinpath(
-    pkgdir(CloudMicrophysics),
-    "src",
-    "parameters",
-    "toml",
-    "SB2006_limiters.toml",
+    pkgdir(CloudMicrophysics), "src", "parameters", "toml", "SB2006_limiters.toml",
 )
 toml_dict = CP.create_toml_dict(FT; override_file)
 SB = CMP.SB2006(toml_dict)
@@ -54,35 +50,35 @@ reff_2M_100_LH  = [CMD.effective_radius_Liu_Hallet_97(cloud_liquid, ρ_air, q_li
 reff_2M_1000_LH = [CMD.effective_radius_Liu_Hallet_97(cloud_liquid, ρ_air, q_liq, FT(1000 * 1e6), FT(0), FT(0)) for q_liq in q_liq_range]
 
 # plotting
-fig = PL.Figure(size = (1100, 1000), fontsize=22, linewidth=3)
+fig = MK.Figure(size = (1100, 1000), fontsize=22, linewidth=3)
 
-ax1 = PL.Axis(fig[1, 1])
-ax2 = PL.Axis(fig[2, 1])
+ax1 = MK.Axis(fig[1, 1])
+ax2 = MK.Axis(fig[2, 1])
 
-PL.ylims!(ax2, [-50, 40])
+MK.ylims!(ax2, [-50, 40])
 
 ax1.xlabel = "q_liq [g/kg]"
 ax1.ylabel = "effective radius [μm]"
 ax2.xlabel = "q_rai [g/kg]"
 ax2.ylabel = "radar reflectivity [dBZ]"
 
-#p_reff_2M_100        = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_100 * 1e6,  color = :blue, linestyle = :dot)
-#p_reff_2M_1000       = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000 * 1e6,  color = :skyblue1, linestyle = :dot)
-p_reff_2M_100_nolim  = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_100_nolim * 1e6,  color = :blue)
-p_reff_2M_100_LH     = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_100_LH * 1e6,     color = :crimson)
+#p_reff_2M_100        = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_100 * 1e6,  color = :blue, linestyle = :dot)
+#p_reff_2M_1000       = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000 * 1e6,  color = :skyblue1, linestyle = :dot)
+p_reff_2M_100_nolim  = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_100_nolim * 1e6,  color = :blue)
+p_reff_2M_100_LH     = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_100_LH * 1e6,     color = :crimson)
 
-p_reff_2M_1000_nolim = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000_nolim * 1e6,  color = :skyblue1)
-p_reff_2M_1000_LH    = PL.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000_LH * 1e6,     color = :orange)
+p_reff_2M_1000_nolim = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000_nolim * 1e6,  color = :skyblue1)
+p_reff_2M_1000_LH    = MK.lines!(ax1, q_liq_range * 1e3,  reff_2M_1000_LH * 1e6,     color = :orange)
 
-p_Z_1M           = PL.lines!(ax2, q_rain_range * 1e3,  Z_1M,  color = :green)
+p_Z_1M           = MK.lines!(ax2, q_rain_range * 1e3,  Z_1M,  color = :green)
 
-p_Z_2M_100       = PL.lines!(ax2, q_rain_range * 1e3,  Z_2M_100,        color = :skyblue1, linestyle = :dot)
-p_Z_2M_100_nolim = PL.lines!(ax2, q_rain_range * 1e3,  Z_2M_100_nolim,  color = :skyblue1)
+p_Z_2M_100       = MK.lines!(ax2, q_rain_range * 1e3,  Z_2M_100,        color = :skyblue1, linestyle = :dot)
+p_Z_2M_100_nolim = MK.lines!(ax2, q_rain_range * 1e3,  Z_2M_100_nolim,  color = :skyblue1)
 
-p_Z_2M_10        = PL.lines!(ax2, q_rain_range * 1e3,  Z_2M_10,        color = :blue, linestyle = :dot)
-p_Z_2M_10_nolim  = PL.lines!(ax2, q_rain_range * 1e3,  Z_2M_10_nolim,  color = :blue)
+p_Z_2M_10        = MK.lines!(ax2, q_rain_range * 1e3,  Z_2M_10,        color = :blue, linestyle = :dot)
+p_Z_2M_10_nolim  = MK.lines!(ax2, q_rain_range * 1e3,  Z_2M_10_nolim,  color = :blue)
 
-PL.Legend(
+MK.Legend(
     fig[1, 2],
     [p_reff_2M_100_nolim, p_reff_2M_100_LH, p_reff_2M_1000_nolim, p_reff_2M_1000_LH],
     [
@@ -93,7 +89,7 @@ PL.Legend(
     ],
     framevisible = false,
 )
-PL.Legend(
+MK.Legend(
     fig[2, 2],
     [p_Z_1M, p_Z_2M_100, p_Z_2M_100_nolim, p_Z_2M_10, p_Z_2M_10_nolim],
     [
@@ -108,4 +104,4 @@ PL.Legend(
 
 #! format: on
 
-PL.save("CloudDiagnostics.svg", fig)
+MK.save("CloudDiagnostics.svg", fig)
diff --git a/docs/src/plots/Frostenberg_fig1.jl b/docs/src/plots/Frostenberg_fig1.jl
index b1aa14baa..9901be893 100644
--- a/docs/src/plots/Frostenberg_fig1.jl
+++ b/docs/src/plots/Frostenberg_fig1.jl
@@ -1,8 +1,6 @@
 import Plots as PL
 
-import CloudMicrophysics as CM
 import CloudMicrophysics.HetIceNucleation as IN
-import ClimaParams as CP
 import CloudMicrophysics.Parameters as CMP
 
 FT = Float32
diff --git a/docs/src/plots/HomFreezingPlots.jl b/docs/src/plots/HomFreezingPlots.jl
index aaac07e47..fbdf763de 100644
--- a/docs/src/plots/HomFreezingPlots.jl
+++ b/docs/src/plots/HomFreezingPlots.jl
@@ -1,17 +1,15 @@
 import CairoMakie as MK
 
 import Thermodynamics as TD
-import CloudMicrophysics as CM
-import ClimaParams as CP
 
-const CMO = CM.Common
-const CMI = CM.HomIceNucleation
-const CMP = CM.Parameters
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.HomIceNucleation as CMI
+import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
-const H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
-const ip = CMP.IceNucleationParameters(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+ip = CMP.IceNucleationParameters(FT)
 
 # Initializing
 T_range = range(229.0, stop = 234.5, length = 50)  # air temperature
@@ -19,15 +17,15 @@ x_sulph = Vector{FT}([0.03, 0.04, 0.06])           # wt% sulphuric acid in dropl
 
 # 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
+    CO.a_w_xT(H2SO4_prs, tps, x_sulph[1], T) - CO.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
+    CO.a_w_xT(H2SO4_prs, tps, x_sulph[2], T) - CO.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
+    CO.a_w_xT(H2SO4_prs, tps, x_sulph[3], T) - CO.a_w_ice(tps, T) for
     T in T_range
 ]
 
@@ -60,7 +58,7 @@ Baum_Delta_a = [0.26, 0.27, 0.28, 0.29, 0.3, 0.32, 0.33, 0.339]
 Baum_J = [4.25e-5, 0.306, 454.09, 2.06e5, 6.31e7, 1.8e12, 4.45e14, 1.69e17]
 
 # Plotting
-fig = MK.Figure(resolution = (800, 500))
+fig = MK.Figure(size = (800, 500))
 ax1 = MK.Axis(
     fig[1, 1],
     ylabel = "log10(J) with J in SI units",
diff --git a/docs/src/plots/Kirkby_organic_nucleation_plots.jl b/docs/src/plots/Kirkby_organic_nucleation_plots.jl
index 1eeace63e..665711ffe 100644
--- a/docs/src/plots/Kirkby_organic_nucleation_plots.jl
+++ b/docs/src/plots/Kirkby_organic_nucleation_plots.jl
@@ -1,8 +1,7 @@
-using Plots
+import Plots
 
-import ClimaParams as CP
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Nucleation as Nucleation
+import CloudMicrophysics.Nucleation
 
 FT = Float64
 params = CMP.OrganicNucleationParameters(FT)
diff --git a/docs/src/plots/KnopfAlpert2013_fig1.jl b/docs/src/plots/KnopfAlpert2013_fig1.jl
index af5b83c65..05173d37c 100644
--- a/docs/src/plots/KnopfAlpert2013_fig1.jl
+++ b/docs/src/plots/KnopfAlpert2013_fig1.jl
@@ -1,15 +1,14 @@
 import Plots as PL
 
 import Thermodynamics as TD
-import CloudMicrophysics as CM
 import CloudMicrophysics.Common as CO
 import CloudMicrophysics.HetIceNucleation as IN
 import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
-const H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
-const kaolinite = CMP.Kaolinite(FT) # dust type
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+kaolinite = CMP.Kaolinite(FT) # dust type
 
 # Initializing
 T_range = range(210, stop = 232, length = 100)  # air temperature
diff --git a/docs/src/plots/KnopfAlpert2013_fig5.jl b/docs/src/plots/KnopfAlpert2013_fig5.jl
index 00d3550db..0bc229996 100644
--- a/docs/src/plots/KnopfAlpert2013_fig5.jl
+++ b/docs/src/plots/KnopfAlpert2013_fig5.jl
@@ -1,9 +1,7 @@
 import CairoMakie as MK
 
-import ClimaParams
 import Thermodynamics as TD
-import CloudMicrophysics as CM
-import CloudMicrophysics.Common as CMO
+import CloudMicrophysics.Common as CO
 import CloudMicrophysics.HetIceNucleation as CMI
 import CloudMicrophysics.Parameters as CMP
 
@@ -40,17 +38,17 @@ T_range = range(228.2, stop = 229.6, length = 100)  # air temperature
 T_dew = FT(228.0)               # dew point temperature
 x_sulph = FT(0)                 # sulphuric acid concentration in droplets
 a_sol = [                       # water activity of solution droplet at equilibrium
-    CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph, T_dew) /
+    CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x_sulph, T_dew) /
     TD.saturation_vapor_pressure(tps, T, TD.Liquid()) for T in T_range
 ]
 # water activity of ice
-a_ice = [CMO.a_w_ice(tps, T) for T in T_range]
+a_ice = [CO.a_w_ice(tps, T) for T in T_range]
 
 Δa_w = @. max(abs(a_sol - a_ice), FT(0.0))
 J_ABIFM = @. CMI.ABIFM_J(illite, Δa_w) * 1e-4 # converted from SI units to cm^-2 s^-1
 
 # Plot results
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 ax1 = MK.Axis(
     fig[1, 1],
     ylabel = "J_het [cm^-2 s^-1]",
diff --git a/docs/src/plots/MarshallPalmer_distribution.jl b/docs/src/plots/MarshallPalmer_distribution.jl
index 9f3c2aaea..4c3ee89d4 100644
--- a/docs/src/plots/MarshallPalmer_distribution.jl
+++ b/docs/src/plots/MarshallPalmer_distribution.jl
@@ -1,14 +1,7 @@
-import CairoMakie as MK
-
-import SpecialFunctions as SF
 
-import Thermodynamics as TD
-import CloudMicrophysics as CM
-import ClimaParams as CP
+import CairoMakie as MK
 
-import Thermodynamics.Parameters as TDP
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Common as CO
 import CloudMicrophysics.Microphysics1M as CM1
 
 FT = Float64
@@ -45,7 +38,7 @@ n_r_2 =
 n_r_3 =
     [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]
 
-fig = MK.Figure(resolution = (1100, 600))
+fig = MK.Figure(size = (1100, 600))
 ax1 = MK.Axis(
     fig[1, 1],
     title = "Marshall-Palmer distribution",
diff --git a/docs/src/plots/Microphysics1M_plots.jl b/docs/src/plots/Microphysics1M_plots.jl
index 73fcaf4cf..eb2cb1704 100644
--- a/docs/src/plots/Microphysics1M_plots.jl
+++ b/docs/src/plots/Microphysics1M_plots.jl
@@ -1,26 +1,21 @@
-import Plots
+import Plots as PL
 
-import Thermodynamics
-import CloudMicrophysics
-import ClimaParams
-import CloudMicrophysics.Common as CO
+import Thermodynamics as TD
 
-const PL = Plots
-const CM1 = CloudMicrophysics.Microphysics1M
-const TD = Thermodynamics
-const CMP = CloudMicrophysics.Parameters
+import CloudMicrophysics.Microphysics1M as CM1
+import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
 
-const tps = Thermodynamics.Parameters.ThermodynamicsParameters(FT)
-const aps = CMP.AirProperties(FT)
-const liquid = CMP.CloudLiquid(FT)
-const ice = CMP.CloudIce(FT)
-const rain = CMP.Rain(FT)
-const snow = CMP.Snow(FT)
-const Chen2022 = CMP.Chen2022VelType(FT)
-const Blk1MVel = CMP.Blk1MVelType(FT)
-const ce = CMP.CollisionEff(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+aps = CMP.AirProperties(FT)
+liquid = CMP.CloudLiquid(FT)
+ice = CMP.CloudIce(FT)
+rain = CMP.Rain(FT)
+snow = CMP.Snow(FT)
+Chen2022 = CMP.Chen2022VelType(FT)
+Blk1MVel = CMP.Blk1MVelType(FT)
+ce = CMP.CollisionEff(FT)
 
 # eq. 5b in [Grabowski1996](@cite)
 function accretion_empirical(q_rai::DT, q_liq::DT, q_tot::DT) where {DT <: Real}
diff --git a/docs/src/plots/Microphysics2M_plots.jl b/docs/src/plots/Microphysics2M_plots.jl
index 1cace8aa7..e85a9d8ee 100644
--- a/docs/src/plots/Microphysics2M_plots.jl
+++ b/docs/src/plots/Microphysics2M_plots.jl
@@ -1,8 +1,8 @@
-using CairoMakie
-CairoMakie.activate!(type = "svg")
+import CairoMakie as MK
+MK.activate!(type = "svg")
 
-import ClimaParams
 import Thermodynamics as TD
+
 import CloudMicrophysics
 import CloudMicrophysics.Microphysics1M as CM1
 import CloudMicrophysics.Microphysics2M as CM2
@@ -10,21 +10,21 @@ import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
 
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
-const aps = CMP.AirProperties(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+aps = CMP.AirProperties(FT)
 
-const KK2000 = CMP.KK2000(FT)
-const B1994 = CMP.B1994(FT)
-const TC1980 = CMP.TC1980(FT)
-const LD2004 = CMP.LD2004(FT)
-const VarTSc = CMP.VarTimescaleAcnv(FT)
-const SB2006 = CMP.SB2006(FT)
+KK2000 = CMP.KK2000(FT)
+B1994 = CMP.B1994(FT)
+TC1980 = CMP.TC1980(FT)
+LD2004 = CMP.LD2004(FT)
+VarTSc = CMP.VarTimescaleAcnv(FT)
+SB2006 = CMP.SB2006(FT)
 
-const ce = CMP.CollisionEff(FT)
+ce = CMP.CollisionEff(FT)
 
-const liquid = CMP.CloudLiquid(FT)
-const rain = CMP.Rain(FT)
-const blk1mvel = CMP.Blk1MVelType(FT)
+liquid = CMP.CloudLiquid(FT)
+rain = CMP.Rain(FT)
+blk1mvel = CMP.Blk1MVelType(FT)
 
 include(
     joinpath(pkgdir(CloudMicrophysics), "docs", "src", "plots", "Wooddata.jl"),
@@ -115,80 +115,59 @@ accK1969_q_rai = [
     q_rai in q_rai_range
 ]
 
-fig = Figure(resolution = (900, 600))
-
-ax1 = Axis(fig[1, 1]; yscale = log10)
-ax2 = Axis(fig[1, 2]; xscale = log10, yscale = log10)
-ax3 = Axis(fig[2, 1]; yscale = log10)
-ax4 = Axis(fig[2, 2]; yscale = log10)
-
-ylims!(ax1, [1e-13, 1e-5])
-ylims!(ax2, [1e-13, 1e-5])
-ylims!(ax3, [1e-7, 5e-6])
-ylims!(ax4, [1e-7, 5e-6])
-
-l1 = lines!(ax1, q_liq_range * 1e3, q_liq_KK2000, color = :red)
-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,
-)
-
-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)
-l12 = lines!(ax2, N_d_range * 1e-6, N_d_TC1980, color = :blue)
-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)
-
-l18 = lines!(ax3, q_liq_range * 1e3, accKK2000_q_liq, color = :red)
-l19 = lines!(ax3, q_liq_range * 1e3, accB1994_q_liq, color = :green)
-l20 = lines!(ax3, q_liq_range * 1e3, accTC1980_q_liq, color = :blue)
-l21 = lines!(ax3, q_liq_range * 1e3, accK1969_q_liq, color = :black)
-
-l22 = lines!(ax4, q_rai_range * 1e3, accKK2000_q_rai, color = :red)
-l23 = lines!(ax4, q_rai_range * 1e3, accB1994_q_rai, color = :green)
-l24 = lines!(ax4, q_rai_range * 1e3, accTC1980_q_rai, color = :blue)
-l25 = lines!(ax4, q_rai_range * 1e3, accK1969_q_rai, color = :black)
-
-l26 = lines!(ax1, q_liq_range * 1e3, q_liq_SB2006, color = :cyan)
-l27 = lines!(ax2, N_d_range * 1e-6, N_d_SB2006, color = :cyan)
-l28 = lines!(ax3, q_liq_range * 1e3, accSB2006_q_liq, color = :cyan)
-l28 = lines!(ax4, q_rai_range * 1e3, accSB2006_q_rai, color = :cyan)
-
-l29 = lines!(ax1, q_liq_range * 1e3, q_liq_VarTimeScaleAcnv, color = :orange)
-l30 = lines!(ax2, N_d_range * 1e-6, N_d_VarTimeScaleAcnv, color = :orange)
-
-ax1.xlabel = "q_liq [g/kg]"
-ax1.ylabel = "autoconversion rate [1/s]"
-ax2.xlabel = "N_d [1/cm3]"
-ax2.ylabel = "autoconversion rate [1/s]"
-ax3.xlabel = "q_liq [g/kg] (q_rai = 0.5 g/kg)"
-ax3.ylabel = "accretion rate [1/s]"
-ax4.xlabel = "q_rai [g/kg] (q_liq = 0.5 g/kg)"
-ax4.ylabel = "accretion rate [1/s]"
-
-Legend(
+fig = MK.Figure(size = (900, 600))
+
+ax1 = MK.Axis(fig[1, 1]; yscale = log10)
+ax2 = MK.Axis(fig[1, 2]; xscale = log10, yscale = log10)
+ax3 = MK.Axis(fig[2, 1]; yscale = log10)
+ax4 = MK.Axis(fig[2, 2]; yscale = log10)
+
+MK.ylims!(ax1, [1e-13, 1e-5]; xlabel = "q_liq [g/kg]", ylabel = "autoconversion rate [1/s]")
+MK.ylims!(ax2, [1e-13, 1e-5]; xlabel = "N_d [1/cm3]", ylabel = "autoconversion rate [1/s]")
+MK.ylims!(ax3, [1e-7, 5e-6]; xlabel = "q_liq [g/kg] (q_rai = 0.5 g/kg)", ylabel = "accretion rate [1/s]")
+MK.ylims!(ax4, [1e-7, 5e-6]; xlabel = "q_rai [g/kg] (q_liq = 0.5 g/kg)", ylabel = "accretion rate [1/s]")
+
+ql_g_kg = q_liq_range * 1e3
+qr_g_kg = q_rai_range * 1e3
+N_d_1_cm3 = N_d_range * 1e-6
+l1 = MK.lines!(ax1, ql_g_kg, q_liq_KK2000, color = :red)
+l2 = MK.lines!(ax1, ql_g_kg, q_liq_B1994, color = :green)
+l3 = MK.lines!(ax1, ql_g_kg, q_liq_TC1980, color = :blue)
+l4 = MK.lines!(ax1, ql_g_kg, q_liq_LD2004, color = :purple)
+l5 = MK.lines!(ax1, ql_g_kg, q_liq_K1969, color = :black)
+l6 = MK.lines!(ax1, KK2000_x_q_liq, KK2000_y_q_liq, color = :red, linestyle = :dash)
+l7 = MK.lines!(ax1, B1994_x_q_liq, B1994_y_q_liq, color = :green, linestyle = :dash)
+l8 = MK.lines!(ax1, TC1980_x_q_liq, TC1980_y_q_liq, color = :blue, linestyle = :dash)
+l9 = MK.lines!(ax1, LD2004_x_q_liq, LD2004_y_q_liq, color = :purple, linestyle = :dash)
+
+l10 = MK.lines!(ax2, N_d_1_cm3, N_d_KK2000, color = :red)
+l11 = MK.lines!(ax2, N_d_1_cm3, N_d_B1994, color = :green)
+l12 = MK.lines!(ax2, N_d_1_cm3, N_d_TC1980, color = :blue)
+l13 = MK.lines!(ax2, N_d_1_cm3, N_d_LD2004, color = :purple)
+l14 = MK.lines!(ax2, KK2000_x_N_d, KK2000_y_N_d, color = :red, linestyle = :dash)
+l15 = MK.lines!(ax2, B1994_x_N_d, B1994_y_N_d, color = :green, linestyle = :dash)
+l16 = MK.lines!(ax2, TC1980_x_N_d, TC1980_y_N_d, color = :blue, linestyle = :dash)
+l17 = MK.lines!(ax2, LD2004_x_N_d, LD2004_y_N_d, color = :purple, linestyle = :dash)
+
+l18 = MK.lines!(ax3, ql_g_kg, accKK2000_q_liq, color = :red)
+l19 = MK.lines!(ax3, ql_g_kg, accB1994_q_liq, color = :green)
+l20 = MK.lines!(ax3, ql_g_kg, accTC1980_q_liq, color = :blue)
+l21 = MK.lines!(ax3, ql_g_kg, accK1969_q_liq, color = :black)
+
+l22 = MK.lines!(ax4, qr_g_kg, accKK2000_q_rai, color = :red)
+l23 = MK.lines!(ax4, qr_g_kg, accB1994_q_rai, color = :green)
+l24 = MK.lines!(ax4, qr_g_kg, accTC1980_q_rai, color = :blue)
+l25 = MK.lines!(ax4, qr_g_kg, accK1969_q_rai, color = :black)
+
+l26 = MK.lines!(ax1, ql_g_kg, q_liq_SB2006, color = :cyan)
+l27 = MK.lines!(ax2, N_d_1_cm3, N_d_SB2006, color = :cyan)
+l28 = MK.lines!(ax3, ql_g_kg, accSB2006_q_liq, color = :cyan)
+l28 = MK.lines!(ax4, qr_g_kg, accSB2006_q_rai, color = :cyan)
+
+l29 = MK.lines!(ax1, ql_g_kg, q_liq_VarTimeScaleAcnv, color = :orange)
+l30 = MK.lines!(ax2, N_d_1_cm3, N_d_VarTimeScaleAcnv, color = :orange)
+
+MK.Legend(
     fig[1, 3],
     [l1, l2, l3, l4, l26, l5, l6, l7, l8, l9, l29],
     [
diff --git a/docs/src/plots/P3AspectRatioPlot.jl b/docs/src/plots/P3AspectRatioPlot.jl
index 24a8808c1..0aef54551 100644
--- a/docs/src/plots/P3AspectRatioPlot.jl
+++ b/docs/src/plots/P3AspectRatioPlot.jl
@@ -1,10 +1,9 @@
-using CairoMakie
-import CloudMicrophysics as CM
-import ClimaParams as CP
+import CairoMakie as MK
+
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.P3Scheme as P3
 
-axis_theme = Theme(
+axis_theme = MK.Theme(
     Axis = (
         xscale = log10,
         xminorticksvisible = true,
@@ -19,7 +18,7 @@ axis_theme = Theme(
     VLines = (linewidth = 1.5,),
 )
 
-logocolors = Makie.Colors.JULIA_LOGO_COLORS
+logocolors = MK.Colors.JULIA_LOGO_COLORS
 cl = [logocolors.blue, logocolors.green, logocolors.red]
 
 #! format: off
@@ -31,20 +30,20 @@ function p3_relations_plot()
 	params = CMP.ParametersP3(FT)
 	
 	D_range = 10.0 .^ (-2:0.01:1.0) * 1e-3
-	fig = Figure(size=(1200, 400), figure_padding = 20)
+	fig = MK.Figure(size=(1200, 400), figure_padding = 20)
 
 	# define plot axis
-	ax1L = Axis(fig[1, 1], ylabel = L"$ϕ(D)$ [-]", title = L"Regimes for $ρ_r = 400$ kg/m³")
-	ax1R = Axis(fig[1, 2], title = L"Regimes for $F_{rim} = 0.4$")
+	ax1L = MK.Axis(fig[1, 1], ylabel = L"$ϕ(D)$ [-]", title = L"Regimes for $ρ_r = 400$ kg/m³")
+	ax1R = MK.Axis(fig[1, 2], title = L"Regimes for $F_{rim} = 0.4$")
 
-	colgap!(fig.layout, 40) # add space between columns so xticklabels don't overlap
+	MK.colgap!(fig.layout, 40) # add space between columns so xticklabels don't overlap
 
 	function lines_and_vlines!(ax, state, color)
-		lines!(ax, D_range * mm, P3.ϕᵢ.(state, D_range); color)
+		MK.lines!(ax, D_range * mm, P3.ϕᵢ.(state, D_range); color)
 
-		vlines!(ax, state.D_th * mm; linestyle = (:dot, :loose), color = :gray)
-		vlines!(ax, state.D_gr * mm; linestyle = :dash, color)
-		vlines!(ax, state.D_cr * mm; linestyle = :dashdot, color)
+		MK.vlines!(ax, state.D_th * mm; linestyle = (:dot, :loose), color = :gray)
+		MK.vlines!(ax, state.D_gr * mm; linestyle = :dash, color)
+		MK.vlines!(ax, state.D_cr * mm; linestyle = :dashdot, color)
 	end
 
 	# Make a plot for each rime fraction
@@ -65,23 +64,23 @@ function p3_relations_plot()
 		lines_and_vlines!(ax1R, state, color)
 	end
 
-	linkyaxes!(ax1R, ax1L)
+	MK.linkyaxes!(ax1R, ax1L)
 
 	# Add legend
 	leg_elems = [
-		map(color -> LineElement(;color), cl), 
-		map(linestyle -> LineElement(;linestyle, linewidth = 1.5, color = :gray), [(:dot, :loose), (:dash, :dense), (:dashdot, :dense)])
+		map(color -> MK.LineElement(;color), cl), 
+		map(linestyle -> MK.LineElement(;linestyle, linewidth = 1.5, color = :gray), [(:dot, :loose), (:dash, :dense), (:dashdot, :dense)])
 	]
 	thresh_labels = [L"$D_{th}$", L"$D_{gr}$", L"$D_{cr}$"]
 	leg_kwargs = (; orientation = :horizontal, nbanks=10, tellheight = false, tellwidth = false, margin = (10, 10, 10, 10), halign = :left, valign = :center)
-	Legend(fig[1,1], leg_elems, [Makie.latexstring.(L"F_{rim} = ", F_rims), thresh_labels], ["Rime fraction", "Thresholds"]; leg_kwargs...)
-	Legend(fig[1,2], leg_elems, [Makie.latexstring.(L"ρ_r = ", ρ_rs), thresh_labels], ["Rime density", "Thresholds"]; leg_kwargs...)
+	MK.Legend(fig[1,1], leg_elems, [MK.latexstring.(L"F_{rim} = ", F_rims), thresh_labels], ["Rime fraction", "Thresholds"]; leg_kwargs...)
+	MK.Legend(fig[1,2], leg_elems, [MK.latexstring.(L"ρ_r = ", ρ_rs), thresh_labels], ["Rime density", "Thresholds"]; leg_kwargs...)
 
-	resize_to_layout!(fig)
+	MK.resize_to_layout!(fig)
 	return fig
 end
 #! format: on
 
-fig = with_theme(p3_relations_plot, axis_theme)
-save("P3Scheme_aspect_ratio.svg", fig)
+fig = MK.with_theme(p3_relations_plot, axis_theme)
+MK.save("P3Scheme_aspect_ratio.svg", fig)
 fig
diff --git a/docs/src/plots/P3ImmersionFreezing.jl b/docs/src/plots/P3ImmersionFreezing.jl
index 9833dfc51..50198d747 100644
--- a/docs/src/plots/P3ImmersionFreezing.jl
+++ b/docs/src/plots/P3ImmersionFreezing.jl
@@ -1,5 +1,5 @@
-import CairoMakie as PL
-PL.activate!(type = "svg")
+import CairoMakie as MK
+MK.activate!(type = "svg")
 
 import Thermodynamics as TD
 import CloudMicrophysics.Parameters as CMP
@@ -61,35 +61,35 @@ dLdt_il_T2 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T2, RH),
 dNdt_il_T2 = [P3.het_ice_nucleation(il, tps, TD.PhasePartition(RH2qₜ(T2, RH), qₗ, qᵢ), Nₗ, RH, T2, p2ρ(T2, RH), dt).dNdt for RH in RH_range]
 
 # plotting
-fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3)
+fig = MK.Figure(size = (1500, 500), fontsize=22, linewidth=3)
 
-ax1 = PL.Axis(fig[1, 1]; yscale = log10)
-ax2 = PL.Axis(fig[1, 2]; yscale = log10)
+ax1 = MK.Axis(fig[1, 1]; yscale = log10)
+ax2 = MK.Axis(fig[1, 2]; yscale = log10)
 
 ax1.xlabel = "RH [%]"
 ax1.ylabel = "ice mass nucleation rate [g/m3/s]"
 ax2.xlabel = "RH [%]"
 ax2.ylabel = "ice number nucleation rate [1/cm3/s]"
 
-l_max_dLdt_T1 = PL.lines!(ax1, RH_range * 1e2,  max_dLdt_T1 * 1e3,  color = :thistle)
-l_max_dLdt_T2 = PL.lines!(ax1, RH_range * 1e2,  max_dLdt_T2 * 1e3,  color = :thistle)
-l_max_dNdt    = PL.lines!(ax2, RH_range * 1e2,  max_dNdt    * 1e-6, color = :thistle)
+l_max_dLdt_T1 = MK.lines!(ax1, RH_range * 1e2,  max_dLdt_T1 * 1e3,  color = :thistle)
+l_max_dLdt_T2 = MK.lines!(ax1, RH_range * 1e2,  max_dLdt_T2 * 1e3,  color = :thistle)
+l_max_dNdt    = MK.lines!(ax2, RH_range * 1e2,  max_dNdt    * 1e-6, color = :thistle)
 
 
-l_dLdt_dd_T1 = PL.lines!(ax1, RH_range * 1e2,  dLdt_dd_T1 * 1e3,  color = :skyblue)
-l_dNdt_dd_T1 = PL.lines!(ax2, RH_range * 1e2,  dNdt_dd_T1 * 1e-6, color = :skyblue)
+l_dLdt_dd_T1 = MK.lines!(ax1, RH_range * 1e2,  dLdt_dd_T1 * 1e3,  color = :skyblue)
+l_dNdt_dd_T1 = MK.lines!(ax2, RH_range * 1e2,  dNdt_dd_T1 * 1e-6, color = :skyblue)
 
-l_dLdt_dd_T2 = PL.lines!(ax1, RH_range * 1e2,  dLdt_dd_T2 * 1e3,  color = :blue3)
-l_dNdt_dd_T2 = PL.lines!(ax2, RH_range * 1e2,  dNdt_dd_T2 * 1e-6, color = :blue3)
+l_dLdt_dd_T2 = MK.lines!(ax1, RH_range * 1e2,  dLdt_dd_T2 * 1e3,  color = :blue3)
+l_dNdt_dd_T2 = MK.lines!(ax2, RH_range * 1e2,  dNdt_dd_T2 * 1e-6, color = :blue3)
 
 
-l_dLdt_il_T1 = PL.lines!(ax1, RH_range * 1e2,  dLdt_il_T1 * 1e3,  color = :orchid)
-l_dNdt_il_T1 = PL.lines!(ax2, RH_range * 1e2,  dNdt_il_T1 * 1e-6, color = :orchid)
+l_dLdt_il_T1 = MK.lines!(ax1, RH_range * 1e2,  dLdt_il_T1 * 1e3,  color = :orchid)
+l_dNdt_il_T1 = MK.lines!(ax2, RH_range * 1e2,  dNdt_il_T1 * 1e-6, color = :orchid)
 
-l_dLdt_il_T2 = PL.lines!(ax1, RH_range * 1e2,  dLdt_il_T2 * 1e3,  color = :purple)
-l_dNdt_il_T2 = PL.lines!(ax2, RH_range * 1e2,  dNdt_il_T2 * 1e-6, color = :purple)
+l_dLdt_il_T2 = MK.lines!(ax1, RH_range * 1e2,  dLdt_il_T2 * 1e3,  color = :purple)
+l_dNdt_il_T2 = MK.lines!(ax2, RH_range * 1e2,  dNdt_il_T2 * 1e-6, color = :purple)
 
-PL.Legend(
+MK.Legend(
     fig[1, 3],
     [l_max_dNdt,
      l_dNdt_dd_T1, l_dNdt_dd_T2,
@@ -103,4 +103,4 @@ PL.Legend(
     ],
     framevisible = false,
 )
-PL.save("P3_het_ice_nucleation.svg", fig)
+MK.save("P3_het_ice_nucleation.svg", fig)
diff --git a/docs/src/plots/P3Melting.jl b/docs/src/plots/P3Melting.jl
index 2f70b929f..8b6e77b5c 100644
--- a/docs/src/plots/P3Melting.jl
+++ b/docs/src/plots/P3Melting.jl
@@ -1,5 +1,5 @@
-import CairoMakie: CairoMakie, Makie
-CairoMakie.activate!(type = "svg")
+import CairoMakie as MK
+MK.activate!(type = "svg")
 
 import Thermodynamics as TD
 import CloudMicrophysics.Parameters as CMP
@@ -56,32 +56,27 @@ dLdt4 = [P3.ice_melt(dist, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt).dL
 dNdt4 = [P3.ice_melt(dist, vel, aps, tps, params.T_freeze .+ ΔT, ρₐ2, dt).dNdt for ΔT in ΔT_range]
 
 # plotting
-fig = Makie.Figure(size = (1500, 500), fontsize=22, linewidth=3)
+fig = MK.Figure(size = (1500, 500), fontsize=22, linewidth=3)
 
-ax1 = Makie.Axis(fig[1, 1]; yscale = log10)
-ax2 = Makie.Axis(fig[1, 2]; yscale = log10)
+ax1 = MK.Axis(fig[1, 1]; yscale = log10, xlabel = "T [C]", ylabel = "ice mass melting rate [g/m³/s]")
+ax2 = MK.Axis(fig[1, 2]; yscale = log10, xlabel = "T [C]", ylabel = "ice number melting rate [1/m³/s]")
 
-ax1.xlabel = "T [C]"
-ax1.ylabel = "ice mass melting rate [g/m3/s]"
-ax2.xlabel = "T [C]"
-ax2.ylabel = "ice number melting rate [1/cm3/s]"
+l_max_dLdt = MK.lines!(ax1, ΔT_range,  max_dLdt * 1e3,  color = :thistle)
+l_max_dNdt = MK.lines!(ax2, ΔT_range,  max_dNdt * 1e-6, color = :thistle)
 
-l_max_dLdt = Makie.lines!(ax1, ΔT_range,  max_dLdt * 1e3,  color = :thistle)
-l_max_dNdt = Makie.lines!(ax2, ΔT_range,  max_dNdt * 1e-6, color = :thistle)
+l_dLdt1 = MK.lines!(ax1, ΔT_range,  dLdt1 * 1e3,  color = :skyblue)
+l_dNdt1 = MK.lines!(ax2, ΔT_range,  dNdt1 * 1e-6, color = :skyblue)
 
-l_dLdt1 = Makie.lines!(ax1, ΔT_range,  dLdt1 * 1e3,  color = :skyblue)
-l_dNdt1 = Makie.lines!(ax2, ΔT_range,  dNdt1 * 1e-6, color = :skyblue)
+l_dLdt2 = MK.lines!(ax1, ΔT_range,  dLdt2 * 1e3,  color = :blue3)
+l_dNdt2 = MK.lines!(ax2, ΔT_range,  dNdt2 * 1e-6, color = :blue3)
 
-l_dLdt2 = Makie.lines!(ax1, ΔT_range,  dLdt2 * 1e3,  color = :blue3)
-l_dNdt2 = Makie.lines!(ax2, ΔT_range,  dNdt2 * 1e-6, color = :blue3)
+l_dLdt3 = MK.lines!(ax1, ΔT_range,  dLdt3 * 1e3,  color = :orchid)
+l_dNdt3 = MK.lines!(ax2, ΔT_range,  dNdt3 * 1e-6, color = :orchid)
 
-l_dLdt3 = Makie.lines!(ax1, ΔT_range,  dLdt3 * 1e3,  color = :orchid)
-l_dNdt3 = Makie.lines!(ax2, ΔT_range,  dNdt3 * 1e-6, color = :orchid)
+l_dLdt4 = MK.lines!(ax1, ΔT_range,  dLdt4 * 1e3,  color = :purple)
+l_dNdt4 = MK.lines!(ax2, ΔT_range,  dNdt4 * 1e-6, color = :purple)
 
-l_dLdt4 = Makie.lines!(ax1, ΔT_range,  dLdt4 * 1e3,  color = :purple)
-l_dNdt4 = Makie.lines!(ax2, ΔT_range,  dNdt4 * 1e-6, color = :purple)
-
-Makie.Legend(
+MK.Legend(
     fig[1, 3],
     [l_max_dNdt, l_dNdt1, l_dNdt2, l_dNdt3, l_dNdt4],
     [
@@ -93,4 +88,4 @@ Makie.Legend(
     ],
     framevisible = false,
 )
-Makie.save("P3_ice_melt.svg", fig)
+MK.save("P3_ice_melt.svg", fig)
diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl
index 3b0c07db3..794a16707 100644
--- a/docs/src/plots/P3SchemePlots.jl
+++ b/docs/src/plots/P3SchemePlots.jl
@@ -1,14 +1,14 @@
-import CairoMakie: Makie
-import CloudMicrophysics as CM
-import ClimaParams as CP
+
+import CairoMakie as MK
+
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.P3Scheme as P3
 
-axis_theme = Makie.Theme(
+axis_theme = MK.Theme(
     Axis = (
         xscale = log10,
         xminorticksvisible = true,
-        xminorticks = Makie.IntervalsBetween(5),
+        xminorticks = MK.IntervalsBetween(5),
         xticks = [0.01, 0.1, 1, 10],
         limits = ((0.01, 10.0), nothing),
         xgridvisible = false,
@@ -18,7 +18,7 @@ axis_theme = Makie.Theme(
     VLines = (linewidth = 1.5,),
 )
 
-logocolors = Makie.Colors.JULIA_LOGO_COLORS
+logocolors = MK.Colors.JULIA_LOGO_COLORS
 cl = [logocolors.blue, logocolors.green, logocolors.red]
 
 #! format: off
@@ -30,29 +30,29 @@ function p3_relations_plot()
 	params = CMP.ParametersP3(FT)
 	
 	D_range = 10.0 .^ (-2:0.01:2.0) * 1e-3
-	fig = Makie.Figure(size=(1200, 900), figure_padding = 20)
+	fig = MK.Figure(size=(1200, 900), figure_padding = 20)
 
 	# define plot axis
-	ax1L = Makie.Axis(fig[1, 1], yscale = log10, ylabel = Makie.L"$m(D)$ (kg)", title = Makie.L"Regimes for $ρ_r = 400$ kg/m³")
-	ax1R = Makie.Axis(fig[1, 2], yscale = log10, title = Makie.L"Regimes for $F_{rim} = 0.95$")
-	ax2L = Makie.Axis(fig[2, 1], yscale = log10, ylabel = Makie.L"$A(D)$ (m²)")
-	ax2R = Makie.Axis(fig[2, 2], yscale = log10)
-	ax3L = Makie.Axis(fig[3, 1], xlabel = Makie.L"$D$ (mm)", ylabel = Makie.L"$ρ$ (kg/m³)")
-	ax3R = Makie.Axis(fig[3, 2], xlabel = Makie.L"$D$ (mm)")
-	map(Makie.contents(fig[:, 2])) do ax ax.yticklabelsvisible = false end
-	map(Makie.contents(fig[1:2, :])) do ax ax.xticklabelsvisible = false end
+	ax1L = MK.Axis(fig[1, 1], yscale = log10, ylabel = MK.L"$m(D)$ (kg)", title = MK.L"Regimes for $ρ_r = 400$ kg/m³")
+	ax1R = MK.Axis(fig[1, 2], yscale = log10, title = MK.L"Regimes for $F_{rim} = 0.95$")
+	ax2L = MK.Axis(fig[2, 1], yscale = log10, ylabel = MK.L"$A(D)$ (m²)")
+	ax2R = MK.Axis(fig[2, 2], yscale = log10)
+	ax3L = MK.Axis(fig[3, 1], xlabel = MK.L"$D$ (mm)", ylabel = MK.L"$ρ$ (kg/m³)")
+	ax3R = MK.Axis(fig[3, 2], xlabel = MK.L"$D$ (mm)")
+	map(MK.contents(fig[:, 2])) do ax ax.yticklabelsvisible = false end
+	map(MK.contents(fig[1:2, :])) do ax ax.xticklabelsvisible = false end
 
-	Makie.colgap!(fig.layout, 40) # add space between columns so xticklabels don't overlap
+	MK.colgap!(fig.layout, 40) # add space between columns so xticklabels don't overlap
 
 	function lines_and_vlines!((ax1, ax2, ax3), state, color)
-		Makie.lines!(ax1, D_range * mm, P3.ice_mass.(state, D_range); color)
-		Makie.lines!(ax2, D_range * mm, P3.ice_area.(state, D_range); color)
-		Makie.lines!(ax3, D_range * mm, P3.ice_density.(state, D_range); color)
+		MK.lines!(ax1, D_range * mm, P3.ice_mass.(state, D_range); color)
+		MK.lines!(ax2, D_range * mm, P3.ice_area.(state, D_range); color)
+		MK.lines!(ax3, D_range * mm, P3.ice_density.(state, D_range); color)
 
 		for ax in [ax1, ax2, ax3]
-			Makie.vlines!(ax, state.D_th * mm; linestyle = (:dot, :loose), color = :gray)
-			Makie.vlines!(ax, state.D_gr * mm; linestyle = :dash, color)
-			Makie.vlines!(ax, state.D_cr * mm; linestyle = :dashdot, color)
+			MK.vlines!(ax, state.D_th * mm; linestyle = (:dot, :loose), color = :gray)
+			MK.vlines!(ax, state.D_gr * mm; linestyle = :dash, color)
+			MK.vlines!(ax, state.D_cr * mm; linestyle = :dashdot, color)
 		end
 	end
 
@@ -74,26 +74,26 @@ function p3_relations_plot()
 		lines_and_vlines!((ax1R, ax2R, ax3R), state, color)
 	end
 
-	Makie.linkyaxes!(ax1R, ax1L)
-	Makie.ylims!(ax1R, 1e-14, 1e-4)
-	Makie.linkyaxes!(ax2R, ax2L)
-	Makie.linkyaxes!(ax3R, ax3L)
+	MK.linkyaxes!(ax1R, ax1L)
+	MK.ylims!(ax1R, 1e-14, 1e-4)
+	MK.linkyaxes!(ax2R, ax2L)
+	MK.linkyaxes!(ax3R, ax3L)
 
 	# Add legend
 	leg_elems = [
-		map(color -> Makie.LineElement(;color), cl), 
-		map(linestyle -> Makie.LineElement(;linestyle, linewidth = 1.5, color = :gray), [(:dot, :loose), (:dash, :dense), (:dashdot, :dense)])
+		map(color -> MK.LineElement(;color), cl), 
+		map(linestyle -> MK.LineElement(;linestyle, linewidth = 1.5, color = :gray), [(:dot, :loose), (:dash, :dense), (:dashdot, :dense)])
 	]
-	thresh_labels = [Makie.L"$D_{th}$", Makie.L"$D_{gr}$", Makie.L"$D_{cr}$"]
+	thresh_labels = [MK.L"$D_{th}$", MK.L"$D_{gr}$", MK.L"$D_{cr}$"]
 	leg_kwargs = (; orientation = :horizontal, nbanks=10, tellheight = false, tellwidth = false, margin = (10, 10, 10, 10), halign = :left, valign = :top)
-	Makie.Legend(fig[1,1], leg_elems, [Makie.latexstring.(Makie.L"F_{rim} = ", F_rims), thresh_labels], ["Rime fraction", "Thresholds"]; leg_kwargs...)
-	Makie.Legend(fig[1,2], leg_elems, [Makie.latexstring.(Makie.L"ρ_r = ", ρ_rs), thresh_labels], ["Rime density", "Thresholds"]; leg_kwargs...)
+	MK.Legend(fig[1,1], leg_elems, [MK.latexstring.(MK.L"F_{rim} = ", F_rims), thresh_labels], ["Rime fraction", "Thresholds"]; leg_kwargs...)
+	MK.Legend(fig[1,2], leg_elems, [MK.latexstring.(MK.L"ρ_r = ", ρ_rs), thresh_labels], ["Rime density", "Thresholds"]; leg_kwargs...)
 
-	Makie.resize_to_layout!(fig)
+	MK.resize_to_layout!(fig)
 	return fig
 end
 #! format: on
 
-fig = Makie.with_theme(p3_relations_plot, axis_theme)
-Makie.save("P3Scheme_relations.svg", fig)
+fig = MK.with_theme(p3_relations_plot, axis_theme)
+MK.save("P3Scheme_relations.svg", fig)
 fig
diff --git a/docs/src/plots/P3SlopeParameterizations.jl b/docs/src/plots/P3SlopeParameterizations.jl
index b91e95af0..e2007f263 100644
--- a/docs/src/plots/P3SlopeParameterizations.jl
+++ b/docs/src/plots/P3SlopeParameterizations.jl
@@ -1,6 +1,6 @@
-import CairoMakie: Makie
-import CloudMicrophysics as CM
-import ClimaParams as CP
+
+import CairoMakie as MK
+
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.P3Scheme as P3
 import Printf
@@ -10,12 +10,12 @@ FT = Float64
 logλs = @. log(10.0^FT(3:0.01:5))
 
 function make_slope_plot(slope_law, title)
-    fig = Makie.Figure(size = (400, 300), figure_padding = 20)
+    fig = MK.Figure(size = (400, 300), figure_padding = 20)
 
-    ax = Makie.Axis(fig[1, 1]; title, xscale = log10, xlabel = "λ", ylabel = "μ")
-    Makie.vlines!(λ_bnds, color = (:gray, 0.5))
+    ax = MK.Axis(fig[1, 1]; title, xscale = log10, xlabel = "λ", ylabel = "μ")
+    MK.vlines!(λ_bnds, color = (:gray, 0.5))
 
-    Makie.lines!(ax, exp.(logλs), P3.get_μ.(slope_law, logλs))
+    MK.lines!(ax, exp.(logλs), P3.get_μ.(slope_law, logλs))
 
     return fig
 end
@@ -28,10 +28,10 @@ slope_power_law = params.slope
 log_λ_from_μ(spl::CMP.SlopePowerLaw, μ) = log((μ + spl.c) / spl.a) / spl.b
 λ_bnds = log_λ_from_μ.(slope_power_law, [0.0, 6.0]) .|> exp
 
-fig = Makie.with_theme(Makie.theme_minimal()) do
+fig = MK.with_theme(MK.theme_minimal()) do
     make_slope_plot(slope_power_law, "μ as a function of λ (power law)")
 end
-Makie.save("P3SlopeParameterizations_power_law.svg", fig)
+MK.save("P3SlopeParameterizations_power_law.svg", fig)
 
 # Multiple solutions issue
 
@@ -46,20 +46,20 @@ function make_multiple_solutions_plot()
     dists = P3.get_distribution_parameters_all_solutions(state; L = L_known, N = N_known)
 
     # Make figure
-    fig = Makie.Figure(size = (800, 300), figure_padding = 20)
-    ax = Makie.Axis(fig[1, 1];
+    fig = MK.Figure(size = (800, 300), figure_padding = 20)
+    ax = MK.Axis(fig[1, 1];
         title = "Valid solutions of shape solver",
         xscale = log10, xlabel = "λ", ylabel = "shape solution", limits = ((10^3.55, 10^4.7), (-0.5, 0.5)),
     )
-    Makie.vlines!(λ_bnds, color = (:gray, 0.5))
-    Makie.hlines!(ax, 0; color = (:gray, 0.2))
+    MK.vlines!(λ_bnds, color = (:gray, 0.5))
+    MK.hlines!(ax, 0; color = (:gray, 0.2))
 
-    Makie.lines!(ax, exp.(logλs), shape_sol; color = :black)
+    MK.lines!(ax, exp.(logλs), shape_sol; color = :black)
     map(dists) do dist
         λ = exp(dist.log_λ)
-        Makie.vlines!(ax, λ; label = Printf.@sprintf("λ = %.0f m⁻¹", λ))
+        MK.vlines!(ax, λ; label = Printf.@sprintf("λ = %.0f m⁻¹", λ))
     end
-    Makie.axislegend(ax; framevisible = true)
+    MK.axislegend(ax; framevisible = true)
 
     # N_ice axis
     mm_to_m = 1e-3 # m
@@ -72,24 +72,24 @@ function make_multiple_solutions_plot()
         xscale = log10,
         xlabel = "D (mm)",
     )
-    Makie.Axis(
+    MK.Axis(
         fig[1, 2];
         ax_opts...,
         # ylabel = L"$N_{ice}$ (1/cm³)", 
         ylabel = "[1/cm⁴]",
         title = "Number concentration distribution, N'(D)",
     )
-    Makie.vlines!(state.D_th * m_to_mm, color = (:gray, 0.5))
+    MK.vlines!(state.D_th * m_to_mm, color = (:gray, 0.5))
     for dist in dists
         N′ = @. exp(P3.log_N′ice(dist, Ds))
-        Makie.lines!(Ds_mm, N′ / m_to_cm^4)
-        # Makie.lines!(Ds_mm, cumsum(N′.(Ds) .* Ds / m_to_cm^3))  # CUMULATIVE DISTRIBUTION
+        MK.lines!(Ds_mm, N′ / m_to_cm^4)
+        # MK.lines!(Ds_mm, cumsum(N′.(Ds) .* Ds / m_to_cm^3))  # CUMULATIVE DISTRIBUTION
     end
 
     return fig
 end
 
-fig = Makie.with_theme(Makie.theme_minimal()) do
+fig = MK.with_theme(MK.theme_minimal()) do
     make_multiple_solutions_plot()
 end
-Makie.save("P3SlopeParameterizations_multiple_solutions.svg", fig)
+MK.save("P3SlopeParameterizations_multiple_solutions.svg", fig)
diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl
index 1d0b831a8..28ed262e0 100644
--- a/docs/src/plots/P3TerminalVelocityPlots.jl
+++ b/docs/src/plots/P3TerminalVelocityPlots.jl
@@ -1,10 +1,8 @@
-import ClimaParams as CP
-import CloudMicrophysics as CM
-import CloudMicrophysics.P3Scheme as P3
-import CloudMicrophysics.Parameters as CMP
-import CairoMakie: Makie
 
-const PSP3 = CMP.ParametersP3
+import CairoMakie as MK
+
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.P3Scheme as P3
 
 FT = Float64
 
@@ -73,7 +71,7 @@ function get_values(
     return (; F_rims, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m, V_m_ϕ)
 end
 
-theme = Makie.Theme(
+theme = MK.Theme(
     Axis = (;
         width = 350,
         height = 350,
@@ -118,34 +116,34 @@ function figure_2()
         get_values(params, Chen2022, L_l, N_l, ρ_a, xres, yres)
 
     ### PLOT ###
-    fig = Makie.Figure()
+    fig = MK.Figure()
 
     # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015
 
-    colormap = Makie.cgrad(:PuBuGn_3, 3, categorical = true)
+    colormap = MK.cgrad(:PuBuGn_3, 3, categorical = true)
     regime_contour_kwargs = (; levels = 3, colormap)
 
     row = 1
-    ax1 = Makie.Axis(fig[row, 1]; title = "Particle regimes with small Dₘ")
-    hm = Makie.contourf!(ax1, F_rims, ρ_rs, D_m_regimes_s; regime_contour_kwargs...)
+    ax1 = MK.Axis(fig[row, 1]; title = "Particle regimes with small Dₘ")
+    hm = MK.contourf!(ax1, F_rims, ρ_rs, D_m_regimes_s; regime_contour_kwargs...)
 
-    ax2 = Makie.Axis(fig[row, 2]; title = "Particle regimes with medium Dₘ")
-    hm = Makie.contourf!(ax2, F_rimm, ρ_rm, D_m_regimes_m; regime_contour_kwargs...)
+    ax2 = MK.Axis(fig[row, 2]; title = "Particle regimes with medium Dₘ")
+    hm = MK.contourf!(ax2, F_rimm, ρ_rm, D_m_regimes_m; regime_contour_kwargs...)
 
-    ax3 = Makie.Axis(fig[row, 3]; title = "Particle regimes with large Dₘ")
-    hm = Makie.contourf!(ax3, F_riml, ρ_rl, D_m_regimes_l; regime_contour_kwargs...)
+    ax3 = MK.Axis(fig[row, 3]; title = "Particle regimes with large Dₘ")
+    hm = MK.contourf!(ax3, F_riml, ρ_rl, D_m_regimes_l; regime_contour_kwargs...)
 
     map(1:3) do col
         ticks = (
             [1, 2, 3],
             ["dense\n nonspherical ice", "graupel", "partially\n rimed ice"],
         )
-        Makie.Colorbar(
+        MK.Colorbar(
             fig[row, col];
             colormap,
             ticks,
             vertical = false,
-            width = Makie.Relative(0.95),
+            width = MK.Relative(0.95),
             height = 10,
             halign = 0.5,
             valign = 0.02,
@@ -165,7 +163,7 @@ function figure_2()
         title = "",
     )
         gp = fig[row, col]
-        ax = Makie.Axis(gp; title)
+        ax = MK.Axis(gp; title)
         row3_opts =
             row == 3 ?
             (;
@@ -176,19 +174,19 @@ function figure_2()
                 bottomspinecolor = :white,
                 topspinecolor = :white,
             ) : (;)
-        hm = Makie.contourf!(ax, F_rim, ρ_r, cfvals)
-        Makie.Colorbar(
+        hm = MK.contourf!(ax, F_rim, ρ_r, cfvals)
+        MK.Colorbar(
             gp,
             hm;
             halign = 0.05,
             valign = 0.05,
-            height = Makie.Relative(0.60),
+            height = MK.Relative(0.60),
             width = 10,
             tellwidth = false,
             ticklabelpad = 0,
             row3_opts...,
         )
-        !isnothing(cvals) && Makie.contour!(ax, F_rim, ρ_r, cvals)
+        !isnothing(cvals) && MK.contour!(ax, F_rim, ρ_r, cvals)
     end
 
     row += 1
@@ -223,14 +221,14 @@ function figure_2()
     title = "Vₘ (using ϕᵢ) with large Dₘ"
     make_plots(row, 3, F_riml, ρ_rl; cfvals = V_m_ϕ_l, cvals = D_m_l, title)
 
-    axs = filter(ax -> ax isa Makie.Axis, fig.content)
-    Makie.linkaxes!(axs...)
+    axs = filter(ax -> ax isa MK.Axis, fig.content)
+    MK.linkaxes!(axs...)
 
-    Makie.resize_to_layout!(fig)
-    Makie.save("MorrisonandMilbrandtFig2.svg", fig)
+    MK.resize_to_layout!(fig)
+    MK.save("MorrisonandMilbrandtFig2.svg", fig)
     fig
 end
 #! format: on
 
 # Terminal Velocity figure
-Makie.with_theme(figure_2, theme)
+MK.with_theme(figure_2, theme)
diff --git a/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl b/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl
index 2a7ebc619..dee08e71e 100644
--- a/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl
+++ b/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl
@@ -1,10 +1,7 @@
-import ClimaParams
-import CloudMicrophysics as CM
+import CairoMakie as MK
+
 import CloudMicrophysics.P3Scheme as P3
 import CloudMicrophysics.Parameters as CMP
-import CairoMakie as Plt
-
-const PSP3 = CMP.ParametersP3
 
 FT = Float64
 
@@ -13,7 +10,7 @@ p3 = CMP.ParametersP3(FT)
 # Testing terminal velocity with liquid fraction
 
 function get_values(
-    p3::PSP3,
+    p3,
     Chen2022::CMP.Chen2022VelType,
     L::FT,
     N::FT,
@@ -53,7 +50,7 @@ function get_values(
 end
 
 function make_axis(fig, row, col, title)
-    return Plt.Axis(
+    return MK.Axis(
         fig[row, col],
         xlabel = "F_rim",
         ylabel = "ρ_r",
@@ -102,7 +99,7 @@ function figure_2()
     (F_rimm_9, ρ_rm_9, V_mm_9, D_mm_9) = get_values(p3, Chen2022, L_m_9, N_m, FT(0.9), ρ_a, xres, yres)
     (F_riml_9, ρ_rl_9, V_ml_9, D_ml_9) = get_values(p3, Chen2022, L_l_9, N_l, FT(0.9), ρ_a, xres, yres)
 
-    fig = Plt.Figure()
+    fig = MK.Figure()
 
     # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015
 
@@ -114,77 +111,77 @@ function figure_2()
     crange_large(len) = range(5, 10, length = len)
 
     ax1 = make_axis(fig, 1, 1, "Vₘ with small Dₘ, F_liq = 0")
-    hm = Plt.contourf!(ax1, F_rims_0, ρ_rs_0, V_ms_0, levels = crange_small(10), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax1, F_rims_0, ρ_rs_0, D_ms_0; args...)
-    # Plt.Colorbar(fig[2, 1], hm, vertical = false)
+    hm = MK.contourf!(ax1, F_rims_0, ρ_rs_0, V_ms_0, levels = crange_small(10), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax1, F_rims_0, ρ_rs_0, D_ms_0; args...)
+    # MK.Colorbar(fig[2, 1], hm, vertical = false)
 
     ax2 = make_axis(fig, 1, 2, "Vₘ with medium Dₘ, F_liq = 0")
-    hm = Plt.contourf!(ax2, F_rimm_0, ρ_rm_0, V_mm_0, levels = crange_med(10), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax2, F_rimm_0, ρ_rm_0, D_mm_0; args...)
-    # Plt.Colorbar(fig[2, 2], hm, vertical = false)
+    hm = MK.contourf!(ax2, F_rimm_0, ρ_rm_0, V_mm_0, levels = crange_med(10), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax2, F_rimm_0, ρ_rm_0, D_mm_0; args...)
+    # MK.Colorbar(fig[2, 2], hm, vertical = false)
 
     ax3 = make_axis(fig, 1, 3, "Vₘ with large Dₘ, F_liq = 0")
-    hm = Plt.contourf!(ax3, F_riml_0, ρ_rl_0, V_ml_0, levels = crange_large(10), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax3, F_riml_0, ρ_rl_0, D_ml_0; args...)
-    # Plt.Colorbar(fig[2, 3], hm, vertical = false)
+    hm = MK.contourf!(ax3, F_riml_0, ρ_rl_0, V_ml_0, levels = crange_large(10), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax3, F_riml_0, ρ_rl_0, D_ml_0; args...)
+    # MK.Colorbar(fig[2, 3], hm, vertical = false)
 
-    Plt.linkxaxes!(ax1, ax2)
-    Plt.linkxaxes!(ax2, ax3)
-    Plt.linkyaxes!(ax1, ax2)
-    Plt.linkyaxes!(ax2, ax3)
+    MK.linkxaxes!(ax1, ax2)
+    MK.linkxaxes!(ax2, ax3)
+    MK.linkyaxes!(ax1, ax2)
+    MK.linkyaxes!(ax2, ax3)
 
     ## Plot F_liq = 0.5 as second row of comparisons
 
     ax4 = make_axis(fig, 3, 1, "Vₘ with small Dₘ, F_liq = 0.5")
-    hm = Plt.contourf!(ax4, F_rims_5, ρ_rs_5, V_ms_5, levels = crange_small(25), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax4, F_rims_5, ρ_rs_5, D_ms_5; args...)
-    # Plt.Colorbar(fig[4, 1], hm, vertical = false)
+    hm = MK.contourf!(ax4, F_rims_5, ρ_rs_5, V_ms_5, levels = crange_small(25), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax4, F_rims_5, ρ_rs_5, D_ms_5; args...)
+    # MK.Colorbar(fig[4, 1], hm, vertical = false)
 
     ax5 = make_axis(fig, 3, 2, "Vₘ with medium Dₘ, F_liq = 0.5")
-    hm = Plt.contourf!(ax5, F_rimm_5, ρ_rm_5, V_mm_5, levels = crange_med(25), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax5, F_rimm_5, ρ_rm_5, D_mm_5; args...)
-    # Plt.Colorbar(fig[4, 2], hm, vertical = false)
+    hm = MK.contourf!(ax5, F_rimm_5, ρ_rm_5, V_mm_5, levels = crange_med(25), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax5, F_rimm_5, ρ_rm_5, D_mm_5; args...)
+    # MK.Colorbar(fig[4, 2], hm, vertical = false)
 
     ax6 = make_axis(fig, 3, 3, "Vₘ with large Dₘ, F_liq = 0.5")
-    hm = Plt.contourf!(ax6, F_riml_5, ρ_rl_5, V_ml_5, levels = crange_large(30), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax6, F_riml_5, ρ_rl_5, D_ml_5; args...)
-    # Plt.Colorbar(fig[4, 3], hm, vertical = false)
+    hm = MK.contourf!(ax6, F_riml_5, ρ_rl_5, V_ml_5, levels = crange_large(30), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax6, F_riml_5, ρ_rl_5, D_ml_5; args...)
+    # MK.Colorbar(fig[4, 3], hm, vertical = false)
 
-    Plt.linkxaxes!(ax1, ax4)
-    Plt.linkxaxes!(ax4, ax5)
-    Plt.linkxaxes!(ax5, ax6)
-    Plt.linkyaxes!(ax1, ax4)
-    Plt.linkyaxes!(ax4, ax5)
-    Plt.linkyaxes!(ax5, ax6)
+    MK.linkxaxes!(ax1, ax4)
+    MK.linkxaxes!(ax4, ax5)
+    MK.linkxaxes!(ax5, ax6)
+    MK.linkyaxes!(ax1, ax4)
+    MK.linkyaxes!(ax4, ax5)
+    MK.linkyaxes!(ax5, ax6)
 
     ## Plot F_liq = 0.9 as second row of comparisons
 
     ax7 = make_axis(fig, 5, 1, "Vₘ with small Dₘ, F_liq = 0.9")
-    hm = Plt.contourf!(ax7, F_rims_9, ρ_rs_9, V_ms_9, levels = crange_small(45), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax7, F_rims_9, ρ_rs_9, D_ms_9; args...)
-    Plt.Colorbar(fig[6, 1], hm, vertical = false)
+    hm = MK.contourf!(ax7, F_rims_9, ρ_rs_9, V_ms_9, levels = crange_small(45), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax7, F_rims_9, ρ_rs_9, D_ms_9; args...)
+    MK.Colorbar(fig[6, 1], hm, vertical = false)
 
     ax8 = make_axis(fig, 5, 2, "Vₘ with medium Dₘ, F_liq = 0.9")
-    hm = Plt.contourf!(ax8, F_rimm_9, ρ_rm_9, V_mm_9, levels = crange_med(45), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax8, F_rimm_9, ρ_rm_9, D_mm_9; args...)
-    Plt.Colorbar(fig[6, 2], hm, vertical = false)
+    hm = MK.contourf!(ax8, F_rimm_9, ρ_rm_9, V_mm_9, levels = crange_med(45), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax8, F_rimm_9, ρ_rm_9, D_mm_9; args...)
+    MK.Colorbar(fig[6, 2], hm, vertical = false)
 
     ax9 = make_axis(fig, 5, 3, "Vₘ with large Dₘ, F_liq = 0.9")
-    hm = Plt.contourf!(ax9, F_riml_9, ρ_rl_9, V_ml_9, levels = crange_large(45), extendlow = :auto, extendhigh = :auto)
-    Plt.contour!(ax9, F_riml_9, ρ_rl_9, D_ml_9; args...)
-    Plt.Colorbar(fig[6, 3], hm, vertical = false)
+    hm = MK.contourf!(ax9, F_riml_9, ρ_rl_9, V_ml_9, levels = crange_large(45), extendlow = :auto, extendhigh = :auto)
+    MK.contour!(ax9, F_riml_9, ρ_rl_9, D_ml_9; args...)
+    MK.Colorbar(fig[6, 3], hm, vertical = false)
 
-    Plt.linkxaxes!(ax1, ax7)
-    Plt.linkxaxes!(ax7, ax8)
-    Plt.linkxaxes!(ax8, ax9)
-    Plt.linkyaxes!(ax1, ax7)
-    Plt.linkyaxes!(ax7, ax8)
-    Plt.linkyaxes!(ax8, ax9)
+    MK.linkxaxes!(ax1, ax7)
+    MK.linkxaxes!(ax7, ax8)
+    MK.linkxaxes!(ax8, ax9)
+    MK.linkyaxes!(ax1, ax7)
+    MK.linkyaxes!(ax7, ax8)
+    MK.linkyaxes!(ax8, ax9)
 
 
-    Plt.resize_to_layout!(fig)
-    # Plt.display(fig)
-    Plt.save("MorrisonandMilbrandtFig2_with_F_liq.svg", fig)
+    MK.resize_to_layout!(fig)
+    # MK.display(fig)
+    MK.save("MorrisonandMilbrandtFig2_with_F_liq.svg", fig)
 end
 
 figure_2()
diff --git a/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
index 55f2d62d8..3c102113d 100644
--- a/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
+++ b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl
@@ -1,10 +1,8 @@
-import ClimaParams
-import CloudMicrophysics as CM
+
+import CairoMakie as MK
+
 import CloudMicrophysics.P3Scheme as P3
 import CloudMicrophysics.Parameters as CMP
-import CairoMakie as Plt
-
-const PSP3 = CMP.ParametersP3
 
 FT = Float64
 
@@ -13,7 +11,7 @@ p3 = CMP.ParametersP3(FT)
 # Testing terminal velocity with liquid fraction
 
 function get_values(
-    p3::PSP3,
+    p3,
     Chen2022::CMP.Chen2022VelType,
     L::FT,
     N::FT,
@@ -75,7 +73,7 @@ function get_values(
 end
 
 function make_axis(fig, row, col, title)
-    return Plt.Axis(
+    return MK.Axis(
         fig[row, col],
         xlabel = "F_rim",
         ylabel = "F_liq",
@@ -126,18 +124,18 @@ function fig1()
     )
 
 
-    fig = Plt.Figure()
+    fig = MK.Figure()
 
     ax1_ = make_axis(fig, 1, 1, "Particle regimes with small Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax1_,
         F_rims,
         F_liqs,
         D_ms_regimes,
         levels = 3,
-        colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true),
+        colormap = MK.cgrad(:PuBuGn_3, 3, categorical = true),
     )
-    Plt.Colorbar(
+    MK.Colorbar(
         fig[2, 1],
         hm,
         vertical = false,
@@ -147,15 +145,15 @@ function fig1()
     )
 
     ax2_ = make_axis(fig, 1, 2, "Particle regimes with medium Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax2_,
         F_rimm,
         F_liqm,
         D_mm_regimes,
         levels = 3,
-        colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true),
+        colormap = MK.cgrad(:PuBuGn_3, 3, categorical = true),
     )
-    Plt.Colorbar(
+    MK.Colorbar(
         fig[2, 2],
         hm,
         vertical = false,
@@ -165,15 +163,15 @@ function fig1()
     )
 
     ax3_ = make_axis(fig, 1, 3, "Particle regimes with large Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax3_,
         F_riml,
         F_liql,
         D_ml_regimes,
         levels = 3,
-        colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true),
+        colormap = MK.cgrad(:PuBuGn_3, 3, categorical = true),
     )
-    Plt.Colorbar(
+    MK.Colorbar(
         fig[2, 3],
         hm,
         vertical = false,
@@ -183,7 +181,7 @@ function fig1()
     )
 
     ax1 = make_axis(fig, 3, 1, "Vₘ with small Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax1,
         F_rims,
         F_liqs,
@@ -192,12 +190,12 @@ function fig1()
         extendlow = :auto,
         extendhigh = :auto,
     )
-    Plt.Colorbar(fig[4, 1], hm, vertical = false)
-    Plt.contour!(ax1, F_rims, F_liqs, D_ms; args...)
+    MK.Colorbar(fig[4, 1], hm, vertical = false)
+    MK.contour!(ax1, F_rims, F_liqs, D_ms; args...)
 
 
     ax2 = make_axis(fig, 3, 2, "Vₘ with medium Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax2,
         F_rimm,
         F_liqm,
@@ -206,12 +204,12 @@ function fig1()
         extendlow = :auto,
         extendhigh = :auto,
     )
-    Plt.Colorbar(fig[4, 2], hm, vertical = false)
-    Plt.contour!(ax2, F_rimm, F_liqm, D_mm; args...)
+    MK.Colorbar(fig[4, 2], hm, vertical = false)
+    MK.contour!(ax2, F_rimm, F_liqm, D_mm; args...)
 
 
     ax3 = make_axis(fig, 3, 3, "Vₘ with large Dₘ")
-    hm = Plt.contourf!(
+    hm = MK.contourf!(
         ax3,
         F_riml,
         F_liql,
@@ -220,15 +218,15 @@ function fig1()
         extendlow = :auto,
         extendhigh = :auto,
     )
-    Plt.Colorbar(fig[4, 3], hm, vertical = false)
-    Plt.contour!(ax3, F_riml, F_liql, D_ml; args...)
+    MK.Colorbar(fig[4, 3], hm, vertical = false)
+    MK.contour!(ax3, F_riml, F_liql, D_ml; args...)
 
 
-    # Plt.linkaxes!(ax1_, ax2_, ax3_)
-    Plt.linkaxes!(ax1, ax2, ax3)
+    # MK.linkaxes!(ax1_, ax2_, ax3_)
+    MK.linkaxes!(ax1, ax2, ax3)
 
-    Plt.resize_to_layout!(fig)
-    Plt.save("P3TerminalVelocity_F_liq_rim.svg", fig)
+    MK.resize_to_layout!(fig)
+    MK.save("P3TerminalVelocity_F_liq_rim.svg", fig)
 end
 
 fig1()
diff --git a/docs/src/plots/P3Thresholds.jl b/docs/src/plots/P3Thresholds.jl
index 95e01cb6f..4d7ec5839 100644
--- a/docs/src/plots/P3Thresholds.jl
+++ b/docs/src/plots/P3Thresholds.jl
@@ -1,6 +1,6 @@
-import CairoMakie: Makie
-import CloudMicrophysics as CM
-import ClimaParams as CP
+
+import CairoMakie as MK
+
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.P3Scheme as P3
 
@@ -15,7 +15,7 @@ get_state2(F_rim, ρ_r) = P3.get_state(params; F_rim, ρ_r)
 
 states = get_state2.(F_rims, ρ_rs')
 
-fig = Makie.Figure(size = (1200, 900), figure_padding = 20)
+fig = MK.Figure(size = (1200, 900), figure_padding = 20)
 
 D_th = round(states[1].D_th * 1e3, digits = 4) # mm
 # Make a plot for each threshold
@@ -26,24 +26,24 @@ for (i, key) in enumerate(thresh_keys)
     else
         ""
     end
-    ax = Makie.Axis(fig[i, 1], xlabel = "F_rim", ylabel = "ρ_r", title = title)
+    ax = MK.Axis(fig[i, 1], xlabel = "F_rim", ylabel = "ρ_r", title = title)
     threshold = getproperty.(states, key) * 1e3
-    hm = Makie.heatmap!(ax, F_rims, ρ_rs, threshold)
-    Makie.Colorbar(fig[i, 2], hm, label = "mm")
+    hm = MK.heatmap!(ax, F_rims, ρ_rs, threshold)
+    MK.Colorbar(fig[i, 2], hm, label = "mm")
 end
 
 # Make a plot for each density
 ρ_keys = (:ρ_g, :ρ_d)
 for (i, key) in enumerate(ρ_keys)
-    ax = Makie.Axis(fig[i, 3], xlabel = "F_rim", ylabel = "ρ_r", title = string(key))
+    ax = MK.Axis(fig[i, 3], xlabel = "F_rim", ylabel = "ρ_r", title = string(key))
     density = if key == :ρ_g
         getproperty.(states, key)
     else # :ρ_d
         P3.get_ρ_d.(params.mass, F_rims, ρ_rs')
     end
-    hm = Makie.heatmap!(ax, F_rims, ρ_rs, density)
-    Makie.Colorbar(fig[i, 4], hm, label = "kg/m³")
+    hm = MK.heatmap!(ax, F_rims, ρ_rs, density)
+    MK.Colorbar(fig[i, 4], hm, label = "kg/m³")
 end
 
-Makie.save("P3Thresholds.svg", fig)
+MK.save("P3Thresholds.svg", fig)
 fig
diff --git a/docs/src/plots/PrecipitationSusceptibilityPlots.jl b/docs/src/plots/PrecipitationSusceptibilityPlots.jl
index a675a1a08..518a4ee87 100644
--- a/docs/src/plots/PrecipitationSusceptibilityPlots.jl
+++ b/docs/src/plots/PrecipitationSusceptibilityPlots.jl
@@ -1,10 +1,11 @@
-using CairoMakie
-CairoMakie.activate!(type = "svg")
+
+import CairoMakie as MK
+MK.activate!(type = "svg")
 
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.PrecipitationSusceptibility as CMPS
 
-const FT = Float64
+FT = Float64
 
 scheme = CMP.SB2006(FT)
 
@@ -37,19 +38,16 @@ acc_rates =
         Ref(N_liq),
     )
 
-fig = Figure()
-
-ax = Axis(fig[1, 1])
+fig = MK.Figure()
 
-ax.xlabel = "q_rai / (q_liq + q_rai)"
-ax.ylabel = "Precipitation susceptibility"
+ax = MK.Axis(fig[1, 1]; xlabel = "q_rai / (q_liq + q_rai)", ylabel = "Precipitation susceptibility")
 
-l1 = lines!(ax, τ, [r.d_ln_pp_d_ln_q_liq for r in aut_rates], color = :red)
-l2 = lines!(ax, τ, [r.d_ln_pp_d_ln_q_rai for r in aut_rates], color = :brown)
-l3 = lines!(ax, τ, [r.d_ln_pp_d_ln_q_liq for r in acc_rates], color = :blue)
-l4 = lines!(ax, τ, [r.d_ln_pp_d_ln_q_rai for r in acc_rates], color = :green)
+l1 = MK.lines!(ax, τ, [r.d_ln_pp_d_ln_q_liq for r in aut_rates], color = :red)
+l2 = MK.lines!(ax, τ, [r.d_ln_pp_d_ln_q_rai for r in aut_rates], color = :brown)
+l3 = MK.lines!(ax, τ, [r.d_ln_pp_d_ln_q_liq for r in acc_rates], color = :blue)
+l4 = MK.lines!(ax, τ, [r.d_ln_pp_d_ln_q_rai for r in acc_rates], color = :green)
 
-Legend(
+MK.Legend(
     fig[1, 2],
     [l1, l2, l3, l4],
     [
@@ -60,4 +58,4 @@ Legend(
     ],
 )
 
-save("Glassmeier-Lohmann_Fig2.svg", fig)
+MK.save("Glassmeier-Lohmann_Fig2.svg", fig)
diff --git a/docs/src/plots/RainEvapoartionSB2006.jl b/docs/src/plots/RainEvapoartionSB2006.jl
index 905d857c1..efecb3e7b 100644
--- a/docs/src/plots/RainEvapoartionSB2006.jl
+++ b/docs/src/plots/RainEvapoartionSB2006.jl
@@ -3,8 +3,6 @@ import CairoMakie as MK
 import SpecialFunctions as SF
 
 import Thermodynamics as TD
-import CloudMicrophysics as CM
-import ClimaParams as CP
 
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.Common as CO
@@ -12,9 +10,9 @@ import CloudMicrophysics.Microphysics2M as CM2
 
 FT = Float64
 
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
-const aps = CMP.AirProperties(FT)
-const SB2006 = CMP.SB2006(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+aps = CMP.AirProperties(FT)
+SB2006 = CMP.SB2006(FT)
 
 function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T)
 
@@ -82,7 +80,7 @@ evap_qᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, q, _qᵣ, ρ, Nᵣ, T).ev
 evap_Nᵣ_3n = [CM2.rain_evaporation(SB2006, aps, tps, q, qᵣ, ρ, _Nᵣ, T).evap_rate_1 for _Nᵣ in Nᵣ_range]
 evap_T_3n = [CM2.rain_evaporation(SB2006, aps, tps, q, qᵣ, ρ, Nᵣ, _T).evap_rate_1 for _T in T_range]
 
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 
 ax1 = MK.Axis(fig[1, 1], xlabel = "q_rain [g/kg]", ylabel = "evap rate [1/cm3/s]")
 ax2 = MK.Axis(fig[2, 1], xlabel = "N_rain [1/cm3]", ylabel = "evap rate [1/cm3/s]")
diff --git a/docs/src/plots/Riccobono_mixed_nucleation_plots.jl b/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
index e97df0b6a..f71c13b56 100644
--- a/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
+++ b/docs/src/plots/Riccobono_mixed_nucleation_plots.jl
@@ -1,7 +1,6 @@
-using Plots
+import Plots as PL
 
-import ClimaParams as CP
-import CloudMicrophysics.Nucleation as Nucleation
+import CloudMicrophysics.Nucleation
 import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
@@ -18,8 +17,8 @@ nucleation_rates = map(bioOxOrg_concentrations) do bioOxOrg_conc
     ) * 1e6
 end
 
-Plots.plot(xaxis = :log, yaxis = :log, lw = 3)
-Plots.plot!(
+PL.plot(xaxis = :log, yaxis = :log, lw = 3)
+PL.plot!(
     bioOxOrg_concentrations,
     nucleation_rates,
     label = "Riccobono parameterization",
@@ -53,7 +52,7 @@ CLOUD_points = [
     (196732207.6712905, 2.987642282932369),
 ]
 
-Plots.plot!(
+PL.plot!(
     CLOUD_points,
     seriestype = :scatter,
     label = "CLOUD data",
@@ -64,4 +63,4 @@ Plots.plot!(
     xticks = [1e6, 1e7, 1e8],
     yticks = [0.01, 0.1, 1, 10],
 )
-Plots.svg("Riccobono_nucleation");
+PL.svg("Riccobono_nucleation");
diff --git a/docs/src/plots/T_vs_wateractivity.jl b/docs/src/plots/T_vs_wateractivity.jl
index c22ecf240..d8e1bb154 100644
--- a/docs/src/plots/T_vs_wateractivity.jl
+++ b/docs/src/plots/T_vs_wateractivity.jl
@@ -1,11 +1,9 @@
 import CairoMakie as MK
 
 import Thermodynamics as TD
-import CloudMicrophysics as CM
-import ClimaParams as CP
 
-const CMO = CM.Common
-const CMP = CM.Parameters
+import CloudMicrophysics.Common as CO
+import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
 tps = TD.Parameters.ThermodynamicsParameters(FT)
@@ -14,8 +12,8 @@ H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
 T_range = range(190, stop = 234, length = 100)
 x = FT(0.1)
 #! format: off
-p_sol_1 = [CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x, T) for T in T_range]     # p_sol for concentration x
-p_sol_0 = [CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T) for T in T_range]   # sat vap pressure over pure liq water using p_sol eqn
+p_sol_1 = [CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, x, T) for T in T_range]     # p_sol for concentration x
+p_sol_0 = [CO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T) for T in T_range]   # sat vap pressure over pure liq water using p_sol eqn
 p_sat_liq = [TD.saturation_vapor_pressure(tps, T, TD.Liquid()) for T in T_range]  # sat vap pressure over pure liq water using TD package
 p_sat_ice = [TD.saturation_vapor_pressure(tps, T, TD.Ice()) for T in T_range]     # sat vap pressure over ice using TD package
 
@@ -26,7 +24,7 @@ a_w_ice_alternate = p_sat_ice ./ p_sol_0   # a_ice if sat vapor pressure over pu
 a_w_ice_μ = [exp((210368 + 131.438*T - (3.32373e6 /T) - 41729.1*log(T))/(8.31441*T)) for T in T_range]  # a_ice using chemical potential parameterization
 
 # Plotting. NOTE: all ".* -1" is only to flip x axis
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 ax1 = MK.Axis(fig[1, 1], title = "Temperature vs Water Activity", ylabel = "T [K]", xlabel = "-a_w", limits = ((-1, -0.4), nothing))
 MK.lines!(ax1, a_w .* -1, T_range, label = "CM default a_w", color = :blue)
 MK.lines!(ax1, a_w_alternate .* -1, T_range, label = "a_w using p(0,T)", linestyle = :dash, color = :blue)
diff --git a/docs/src/plots/TerminalVelocity.jl b/docs/src/plots/TerminalVelocity.jl
index 4fa6caecb..ecd871737 100644
--- a/docs/src/plots/TerminalVelocity.jl
+++ b/docs/src/plots/TerminalVelocity.jl
@@ -1,30 +1,28 @@
 import Plots as PL
-using Measures
+import Plots.Measures.mm
 
-import ClimaParams as CP
-
-FT = Float64
 
 import CloudMicrophysics.MicrophysicsNonEq as CMNe
 import CloudMicrophysics.Microphysics1M as CM1
-import CloudMicrophysics.Microphysics2M as CM2
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Common as CMO
+import CloudMicrophysics.Common as CO
+
+FT = Float64
 
-const rain = CMP.Rain(FT)
-const liquid = CMP.CloudLiquid(FT)
-const ice = CMP.CloudIce(FT)
-const snow = CMP.Snow(FT)
+rain = CMP.Rain(FT)
+liquid = CMP.CloudLiquid(FT)
+ice = CMP.CloudIce(FT)
+snow = CMP.Snow(FT)
 
-const SB2006 = CMP.SB2006(FT)
-const SB2006_no_lim = CMP.SB2006(FT, false)
+SB2006 = CMP.SB2006(FT)
+SB2006_no_lim = CMP.SB2006(FT, false)
 
-const Chen2022 = CMP.Chen2022VelType(FT)
-const SB2006Vel = CMP.SB2006VelType(FT)
-const Blk1MVel = CMP.Blk1MVelType(FT)
+Chen2022 = CMP.Chen2022VelType(FT)
+SB2006Vel = CMP.SB2006VelType(FT)
+Blk1MVel = CMP.Blk1MVelType(FT)
 
-const oblate = CM1.Oblate()
-const prolate = CM1.Prolate()
+oblate = CM1.Oblate()
+prolate = CM1.Prolate()
 
 function aspect_ratio_snow_1M_oblate(snow::CMP.Snow, D::FT) where {FT <: Real}
     (; r0, m0, me, χm, Δm) = snow.mass
@@ -62,7 +60,7 @@ function rain_terminal_velocity_individual_Chen(
     ρₐ::FT,
     D::FT, #in m
 ) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B1(velo_scheme, ρₐ)
+    ai, bi, ci = CO.Chen2022_vel_coeffs_B1(velo_scheme, ρₐ)
 
     v = 0
     for i in 1:3
@@ -77,7 +75,7 @@ function ice_terminal_velocity_individual_Chen(
     ρₐ::FT,
     D::FT, #in m
 ) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B2(velo_scheme, ρₐ, ice.ρᵢ)
+    ai, bi, ci = CO.Chen2022_vel_coeffs_B2(velo_scheme, ρₐ, ice.ρᵢ)
 
     v = 0
     for i in 1:2
@@ -92,7 +90,7 @@ function snow_terminal_velocity_individual_Chen(
     ρₐ::FT,
     D::FT, #in m
 ) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
+    ai, bi, ci = CO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
     (; ϕ, κ) = snow.aspr
 
     v = 0
@@ -108,7 +106,7 @@ function snow_terminal_velocity_individual_Chen_oblate(
     ρₐ::FT,
     D_r::FT, #in m
 ) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
+    ai, bi, ci = CO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
 
     ϕ = aspect_ratio_snow_1M_oblate(snow, D_r)
 
@@ -125,7 +123,7 @@ function snow_terminal_velocity_individual_Chen_prolate(
     ρₐ::FT,
     D_r::FT, #in m
 ) where {FT <: Real}
-    ai, bi, ci = CMO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
+    ai, bi, ci = CO.Chen2022_vel_coeffs_B4(velo_scheme, ρₐ, snow.ρᵢ)
 
     ϕ = aspect_ratio_snow_1M_prolate(snow, D_r)
 
diff --git a/docs/src/plots/TerminalVelocity2M.jl b/docs/src/plots/TerminalVelocity2M.jl
index c2cc6b602..234975a70 100644
--- a/docs/src/plots/TerminalVelocity2M.jl
+++ b/docs/src/plots/TerminalVelocity2M.jl
@@ -1,7 +1,6 @@
 import Plots as PL
-using Measures
+import Plots.Measures.mm
 
-import ClimaParams as CP
 import CloudMicrophysics.Microphysics2M as CM2
 import CloudMicrophysics.Parameters as CMP
 
diff --git a/docs/src/plots/Thersholds_transitions.jl b/docs/src/plots/Thersholds_transitions.jl
index 0992a78fe..18478eb01 100644
--- a/docs/src/plots/Thersholds_transitions.jl
+++ b/docs/src/plots/Thersholds_transitions.jl
@@ -1,12 +1,10 @@
-import Plots
-import CloudMicrophysics
-import ClimaParams
 
-const PL = Plots
-const CM1 = CloudMicrophysics.Microphysics1M
-const CM2 = CloudMicrophysics.Microphysics2M
-const CP = ClimaParams
-const CMP = CloudMicrophysics.Parameters
+import Plots as PL
+
+import ClimaParams as CP
+import CloudMicrophysics.Microphysics1M as CM1
+import CloudMicrophysics.Microphysics2M as CM2
+import CloudMicrophysics.Parameters as CMP
 
 FT = Float64
 
diff --git a/docs/src/plots/activity_based_deposition.jl b/docs/src/plots/activity_based_deposition.jl
index 551e16467..a42180425 100644
--- a/docs/src/plots/activity_based_deposition.jl
+++ b/docs/src/plots/activity_based_deposition.jl
@@ -1,16 +1,15 @@
 import Plots as PL
 
 import Thermodynamics as TD
-import CloudMicrophysics as CM
-import CloudMicrophysics.Common as CO
+
 import CloudMicrophysics.HetIceNucleation as IN
 import CloudMicrophysics.Parameters as CMP
 
 FT = Float32
-const tps = TD.Parameters.ThermodynamicsParameters(FT)
-const feldspar = CMP.Feldspar(FT)
-const ferrihydrite = CMP.Ferrihydrite(FT)
-const kaolinite = CMP.Kaolinite(FT)
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+feldspar = CMP.Feldspar(FT)
+ferrihydrite = CMP.Ferrihydrite(FT)
+kaolinite = CMP.Kaolinite(FT)
 
 # Initializing
 Δa_w = range(FT(0), stop = FT(0.32), length = 50)    # difference in solution and ice water activity
diff --git a/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl b/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
index eb14f5dde..04bc81698 100644
--- a/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
+++ b/docs/src/plots/compare_vehkamaki_CLOUD_nucleation.jl
@@ -1,9 +1,7 @@
-using Plots
+import Plots
 
-import ClimaParams as CP
-import CloudMicrophysics as CM
 import CloudMicrophysics.Parameters as CMP
-import CloudMicrophysics.Nucleation as Nucleation
+import CloudMicrophysics.Nucleation
 
 FT = Float64
 params = CMP.H2S04NucleationParameters(FT)
diff --git a/docs/src/plots/lehtinen_apparent_nucleation_rate.jl b/docs/src/plots/lehtinen_apparent_nucleation_rate.jl
index 597b9d098..ce96bbb1a 100644
--- a/docs/src/plots/lehtinen_apparent_nucleation_rate.jl
+++ b/docs/src/plots/lehtinen_apparent_nucleation_rate.jl
@@ -1,4 +1,4 @@
-using Plots
+import Plots
 
 import CloudMicrophysics.Nucleation as Nucleation
 
@@ -20,7 +20,7 @@ for (cs, od) in zip(coag_sinks, output_diams)
     append!(rates, rate)
 end
 
-plot(
+Plots.plot(
     collect(output_diams),
     rates,
     # yaxis = :log,
@@ -28,4 +28,4 @@ plot(
     ylabel = "Apparent Nucleation Rate Reduction",
     legend = false,
 )
-savefig("apparent_nucleation.svg");
+Plots.savefig("apparent_nucleation.svg");
diff --git a/docs/src/plots/linear_HOM_J.jl b/docs/src/plots/linear_HOM_J.jl
index 997c272c2..ac747679a 100644
--- a/docs/src/plots/linear_HOM_J.jl
+++ b/docs/src/plots/linear_HOM_J.jl
@@ -1,10 +1,10 @@
 import CairoMakie as MK
-import CloudMicrophysics as CM
+
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.HomIceNucleation as CMH
 
 FT = Float32
-const ip = CMP.IceNucleationParameters(FT)
+ip = CMP.IceNucleationParameters(FT)
 
 # Koop2000 parameterization
 Koop_Δa = collect(0.26:0.0025:0.34)
@@ -18,7 +18,7 @@ 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)
+fig = MK.Figure(size = (800, 600), fontsize = 18)
 ax1 = MK.Axis(
     fig[1, 1],
     xlabel = "Δa_w [-]",
diff --git a/parcel/Example_P3_vs_activitybased.jl b/parcel/Example_P3_vs_activitybased.jl
index 240eddf48..7003caa20 100644
--- a/parcel/Example_P3_vs_activitybased.jl
+++ b/parcel/Example_P3_vs_activitybased.jl
@@ -69,7 +69,7 @@ deposition_growth = "Deposition"
 size_distribution = "Monodisperse"
 
 # Plotting
-fig = MK.Figure(resolution = (1000, 600))
+fig = MK.Figure(size = (1000, 600))
 ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
 ax2 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]", yscale = log10)
 ax3 = MK.Axis(fig[1, 3], ylabel = "Temperature [K]")