diff --git a/Project.toml b/Project.toml
index 677f5ae81d..885f9808b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,8 +7,8 @@ version = "0.11.1"
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
@@ -17,7 +17,7 @@ CLIMAParameters = "0.7"
 CUDAKernels = "0.4"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"
-KernelAbstractions = "0.8"
+KernelAbstractions = "0.8, 0.9"
 SpecialFunctions = "1, 2"
 Thermodynamics = "0.9, 0.10"
 julia = "1.5"
diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 57f96124cd..a7dca46de3 100644
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -216,6 +216,19 @@ @article{Kirkby2016
   year = {2016}
 }
 
+@Article{KnopfAlpert2013,
+  author = "Knopf, Daniel A. and Alpert, Peter A.",
+  title  = "A water activity based model of heterogeneous ice nucleation kinetics for freezing of water and aqueous solution droplets",
+  journal = "Faraday Discuss.",
+  year = "2013",
+  volume = "165",
+  issue = "0",
+  pages = "513-534",
+  publisher = "The Royal Society of Chemistry",
+  doi = "10.1039/C3FD00035D",
+  url = "http://dx.doi.org/10.1039/C3FD00035D",
+}
+
 @article{Koop2000,
 author = {Koop, Thomas and al, et},
 title = {Water activity as the determinant for homogeneous ice nucleation in aqueous solutions},
@@ -439,14 +452,15 @@ @article{TripoliCotton1980
   doi = {10.1175/1520-0450(1980)019<1037:ANIOSF>2.0.CO;2}
 }
 
-@article{Tully2022,
+@article{Tully2023,
+  title = {Assessing predicted cirrus ice properties between two deterministic ice formation parameterizations},
   author = {Tully, C. and Neubauer, D. and Lohmann, U.},
-  title = {Technical Note: assessing predicted cirrus ice properties between two deterministic ice formation parameterizations},
-  journal = {EGUsphere},
-  volume = {2022},
-  year = {2022},
-  pages = {1--25},
-  doi = {10.5194/egusphere-2022-1057}
+  journal = {Geoscientific Model Development},
+  volume = {16},
+  year = {2023},
+  number = {10},
+  pages = {2957--2973},
+  doi = {10.5194/gmd-16-2957-2023}
 }
 
 @article{Vehkamaki2002,
@@ -492,3 +506,15 @@ @article{SeifertBeheng2006
   year={2006},
   publisher={Springer}
 }
+
+@Article{Spichtinger2023,
+  AUTHOR = {Spichtinger, P. and Marschalik, P. and Baumgartner, M.},
+  TITLE = {Impact of formulations of the homogeneous nucleation rate on ice nucleation events in cirrus},
+  JOURNAL = {Atmospheric Chemistry and Physics},
+  VOLUME = {23},
+  YEAR = {2023},
+  NUMBER = {3},
+  PAGES = {2035--2060},
+  URL = {https://acp.copernicus.org/articles/23/2035/2023/},
+  DOI = {10.5194/acp-23-2035-2023}
+}
diff --git a/docs/make.jl b/docs/make.jl
index 829f3832ca..1580257f10 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -21,6 +21,7 @@ pages = Any[
     "Precipitation Susceptibility" => "PrecipitationSusceptibility.md",
     "API" => "API.md",
     "References" => "References.md",
+    "Developer's Guide" => "DevelopersGuide.md",
 ]
 
 mathengine = MathJax(
diff --git a/docs/src/DevelopersGuide.md b/docs/src/DevelopersGuide.md
new file mode 100644
index 0000000000..ead7f723ca
--- /dev/null
+++ b/docs/src/DevelopersGuide.md
@@ -0,0 +1,194 @@
+# Developer's Guide
+
+This doc is meant to guide any beginner developers with organization
+  and common surprises associated with contributing to `CloudMicrophysics.jl`.
+
+## GitHub workflow
+
+We work on branches of the main repository, rather than personal forks.
+Each branch name should ideally start with your initials, followed by short
+  name relevant to what will be added in that branch.
+It is considered a good practice to open an issue that describes the planned work
+  before starting the implementation.
+This allows you to get feedback from other developers and advertise the planned work to the group.
+When the implemented code is ready for review, create a pull request (PR) on GitHub
+  and tag the issue it is addressing.
+In the PR you can describe which parts of the relevant issue are solved and what needs to
+  be done to reach that goal.
+After creating a PR on GitHub, you will find that every following commit & push you make
+  will go through continuous integration (CI) checks:
+
+- The, `ci / ci 1.8.1 ` runs the tests on Ubuntu, Windows and OSX machines
+  using GitHub Actions.
+  The tests include some unit tests and very simple performance tests.
+  It may happen that the tests pass on your local machine,
+  but fail on one of the machines provided in the cloud for the CI.
+  This is especially true for the performance tests,
+  as the execution times may vary a lot depending on the machine.
+  The CI is the source of truth over local tests, and all tests must pass in the CI before code can be merged.
+  See the [Tests](https://clima.github.io/CloudMicrophysics.jl/dev/DevelopersGuide/#Tests)
+  section below for debugging tips.
+- `CloudMicrophysics.jl` has a small set of tests that are run on the GPUs using `buildkite`.
+  They are triggered automatically when trying to merge a PR,
+  or manually by commenting `bors try` on your PR.
+  If the GPU tests fail, check first for use of any "out of place" functions
+  (for example, `@warn` will not work on GPU).
+- The `Documentation / docbuild` builds the documentation.
+  Most common errors in the docbuild are related to missing docstrings,
+  see the [Documentation](https://clima.github.io/CloudMicrophysics.jl/dev/DevelopersGuide/#Documentation) section for hints.
+  This check may fail if the source code itself is not compiling,
+  so please handle the compilation errors first.
+  If the documentation build was successful, the `documenter/deploy`
+  will display the documentation page based on the PR (click on the details).
+- The `JuliaFormatter / format` ensures consistent formatting throughout the repository.
+  If this check fails, you can click on details
+  to see which files are not following the formatting rules.
+  You can apply the formatter by running
+  `julia --project=.dev .dev/climaformat.jl file_name` in the terminal.
+  The formatter might break if the code does not compile,
+  so make sure you handle the compilation errors first.
+- The `codecov` shows what percentage of source code lines are exercised inside the tests.
+  We strive to keep the test coverage high, but this test is not strictly required to pass before merging a PR.
+
+Once the PR is opened, you can request reviewers to look over and approve your work.
+To keep things tidy, we want PRs to have very few (ideally only one) commit
+  before merging to the main branch.
+You can squash and rebase multiple commits into one
+  in your favorite editor (VS Code, Vim etc).
+In case of doubt, see Git tutorials on how to squash and rebase
+  or reach out to other developers in our team.
+The first couple of times it pays off to create a "backup branch" before starting your rebase,
+  and then comparing afterwards if the rebased and backup branches are identical.
+If the main branch has been updated after your branch was created, 
+you will need to rebase onto the main branch. To do this, 
+run git rebase origin/main and solve any conflicts manually. 
+This is done more easily if you only have one commit.
+
+We (still) use `bors` bot to manage PR queues.
+To merge your PR into the `main` branch, comment `bors r+` and wait for all checks to pass.
+
+The easiest way to use the new additions in the `main` branch of `CloudMicrophysics.jl`
+  from another package is to do a package release.
+To do that, you need to change the package version in the `Project.toml`.
+We do a patch release if the API did not change (for example `0.11.1 -> 0.11.2`)
+  and a minor version release if it did (for example `0.11.1 -> 0.12.0`).
+We use the `JuliaRegistrator` bot to register new versions in the `Julia` package ecosystem,
+  by commenting `@JuliaRegistrator register` under the merged PR that changes the package version.
+You can also use specific branches from the repository,
+  if you are working with commits that were not yet merged into `main`.
+See [Pkg.jl docs](https://pkgdocs.julialang.org/v1/managing-packages/#Adding-registered-packages)
+  for more details.
+
+## New contributions
+
+`CloudMicrophysics.jl` is a collection of point-wise functions grouped
+  in different modules depending on which aspect of aerosol and cloud microphysics
+  they address.
+When adding a new function, decide if a new module is required, or add to the existing one.
+  If the added function will be used in other modules, export the function. Avoid exporting
+  functions that are only used within the module it is defined in or functions that have the
+  same name as a function already in the API.
+Avoid shortening words or phrases when naming new functions. The name of the function
+  should be self-explanatory yet brief.
+All functions should have docstrings describing the API, as well as
+  documentation focusing on the scientific aspects of what they do.
+All functions should have their own unit, performance and GPU tests.
+All free parameters should be stored in [CLIMAParameters](https://github.com/CliMA/CLIMAParameters.jl).
+It is usually faster to prototype defining the free parameters locally,
+  and move them to `CLIMAParameters` at a last step.
+
+Other files that may require editing after you make a new function are:
+ - `CloudMicrophysics.jl` (found in `src` folder) if you need to include a new source file,
+ - `index.md` (found in `docs/src` folder) if you want to mention the new additions in the main documentation page,
+ - `make.jl` (found in `docs` folder) if you are adding a new file to the documentation,
+ - `API.md` (found in `docs/src` folder) if you are adding functions (see ``Documentation`` section for more details),
+ - `runtests.jl` (found in `test` folder) if you are adding new test file for unit tests.
+
+## Documentation
+
+Each new addition to the library should be accompanied by a documentation page
+  summarizing the derivation/assumptions and it's potential uses.
+If possible, it's really appreciated to also add short code snippets that reproduce
+  results from the literature.
+Those code snippets are executed every time documentation is build and provide
+  great examples on how different available parameterizations work.
+
+Additionally, each function in the source code is required to have a docstring
+  and should be added to the `API` documentation page.
+The docstring formatting is pretty strict:
+  (i) no empty lines between the docstring and the function/module,
+  (ii) docstring starts and ends with a line consisting of three quotation marks,
+  (iii) the second line has the function's name preceded by exactly 4 spaces.
+Examples can be found in the source files.
+Be sure to add the function to `API.md` in the docs folder.
+Missing docstrings and functions in the API cause the documentation build to break.
+
+To build and work the documentation locally, you can use [LiveServer.jl](https://github.com/tlienart/LiveServer.jl#serve-docs). This will compile the documentation to a local server that is updated whenever you make changes.
+Alternatively, you can save the documentation to a static webpage: `julia --project=docs docs/make.jl`.
+The index page will be saved in `docs/build` folder.
+
+## Tests
+
+### Unit Tests
+
+Unit tests aim to ensure that parameterizations make physical sense.
+They can be found in the `test` folder under files named after corresponding source files.
+If you create a new function, please also create a new test that checks it.
+If creating a new file for unit tests, make sure you import `Test` and any other necessary libraries.
+There is some boilerplate code needed to create the set with free parameters
+  based on the default `toml` file from [CLIMAParameters](https://github.com/CliMA/CLIMAParameters.jl).
+  This would be the `include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))` line found in
+  the existing unit tests.
+
+Some possible tests include checking if the returned values agree with values
+  in the literature, if something is smaller/greater at warmer/cooler
+  temperatures, if assertion errors are returned when a function is used outside its
+  valid range of parameters, or if a function is zero at certain input values.
+In general, writing good tests is difficult and we are always on the lookout for new good candidates.
+We strive to exercise all functions in some way in tests,
+  so that at minimum we can catch changes in the API.
+
+You can run the tests locally: `julia --project=test test/runtests.jl`.
+
+### Performance Tests
+
+Performance tests check the memory allocations (there should be none) and execution times
+  of some of the functions.
+They are found in the `test` folder under a single file named `performance_tests.jl`.
+
+You can add performance tests for your new functions in the `benchmark_test(FT)` function.
+Parameters are listed at the start of the function.
+Next, add your performance tests under the comment with the file name which your
+  new function is in.
+This is done by calling `bench_press(function_name, (Parameter1, Parameter2), min_time)`.
+The last argument is an estimate of the minimum run-time of the function.
+It may take some trial and error to find a number
+  that will satisfy the `ci` tests that are run on different operating systems
+  and for different floating-point precision.
+You can identify which specific performance test is failing on GitHub
+  by clicking on details next to a failed check from one of the `ci / ci 1.8.1` checks.
+Adjust the estimated minimum run-time appropriately.
+
+### GPU Tests
+
+Graphic Processing Unit (GPU) tests check that `CloudMicrophysics.jl` functions are able to run on GPU.
+They are as simple as checking that a certain input returns a known value.
+Right now, we do not test the whole library on the GPUs,
+  so the support is limited.
+GPU tests can be found in the `test` folder under a single file named `gpu_tests.jl`.
+There are two things to add: a kernel for your function and the actual test.
+
+Kernels are added at the top half of the file using `@kernel`.
+The naming convention follows `test_<function name>_kernel!`.
+Within the kernel, use `@inbounds` to place any outputs into an output array.
+
+The test itself should be added in the `test_gpu(FT)` function.
+The test starts with defining `data_length` and ends after an `@test` macro.
+`data_length` corresponds to the number of outputs your function has.
+Add an array for each input required by your function.
+For example, if you want to test at temperature of 230K,
+  you can add `T = ArrayType([FT(230)])`.
+Add a comment of what you are testing and use `@test` to create your test.
+The GPU tests are ran twice: for `Float64` and `Float32`.
+Similar as with performance tests, some trial and error is needed
+  to find good tolerances for both options.
diff --git a/docs/src/HomFreezingPlots.jl b/docs/src/HomFreezingPlots.jl
new file mode 100644
index 0000000000..c21c3c9804
--- /dev/null
+++ b/docs/src/HomFreezingPlots.jl
@@ -0,0 +1,85 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+const CMT = CM.CommonTypes
+const CMO = CM.Common
+const CMI_het = CM.HetIceNucleation
+const CMI_hom = CM.HomIceNucleation
+const CMP = CM.Parameters
+
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+FT = Float64
+toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+const prs = cloud_microphysics_parameters(toml_dict)
+thermo_params = CMP.thermodynamics_params(prs)
+
+# Initializing
+temp = collect(230.0:0.25:239.5)   # air temperature
+x_sulph = 0.027                    # wt% sulphuric acid in droplets
+Delta_a = Vector{Float64}(undef, length(temp))
+J = Vector{Float64}(undef, length(temp))
+
+it = 1
+for T in temp
+    Delta_a[it] = CMO.Delta_a_w(prs, x_sulph, T)
+    J[it] = CMI_hom.homogeneous_J(Delta_a[it])
+    global it += 1
+end
+log10J = @. log10(J)
+
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(
+    fig[1, 1],
+    ylabel = "log10(J) with J in SI units",
+    xlabel = "Temperature (K)",
+    title = "Homogeneous Ice Nucleation Rate",
+)
+
+Koop2000_temp = [
+    230.0139,
+    231.0179,
+    231.9860,
+    233.00797,
+    234.01195,
+    235.01594,
+    236.00199,
+    237.16733,
+    238.51195,
+    239.51594,
+    240.46613,
+    241.30876,
+    241.91832,
+    242.33068,
+    242.74303,
+    243.10159,
+]
+Koop2000_log10J = [
+    24.33735,
+    22.6506,
+    21.0843,
+    19.5181,
+    18.072289,
+    16.626506,
+    15.24096,
+    13.493975,
+    11.26506,
+    9.39759,
+    7.349397,
+    5.3012,
+    3.674698,
+    2.4698,
+    1.1445783,
+    0,
+]
+
+MK.lines!(ax1, temp, log10J, label = "CliMA")
+MK.lines!(ax1, Koop2000_temp, Koop2000_log10J, label = "Koop 2000")
+
+# MK.Legend(fig[1, 1], ["CliMA", "Koop 2000"])
+fig[1, 2] = MK.Legend(fig, ax1, "Legend", framevisible = false)
+
+MK.save("homogeneous_J.svg", fig)
diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
index 9d34a484ba..30eb7c6ab1 100644
--- a/docs/src/IceNucleation.md
+++ b/docs/src/IceNucleation.md
@@ -115,7 +115,9 @@ The nucleation rate coefficient is determined with the cubic function from [Koop
 ```
 This parameterization is valid only when ``0.26 < \Delta a_w < 0.36`` and ``185K < T < 235K``.
 
-## ABIFM Example Figures
+## Example Figures
+### ABIFM Ice Nucleation Rate Coefficient
+Plotted below is the ABIFM derived ice nucleation rate alongside data and a parameterization from [KnopfAlpert2013](@cite).
 ```@example
 import Plots
 
@@ -170,3 +172,17 @@ PL.plot!(KA13_Delta_a_param, KA13_log10J_param, linecolor = :red, label="paper p
 PL.savefig("Knopf_Alpert_fig_1.svg")
 ```
 ![](Knopf_Alpert_fig_1.svg)
+
+### Homogeneous Ice Nucleation Rate Coefficient
+Here is a comparison of our parameterization of ``J_{hom}`` compared to Koop 2000 as
+  plotted in figure 1 of [Spichtinger2023](@cite). Our parameterization differs in the calculation
+  of ``\Delta a_w``. We define water activity to be a ratio of saturated vapor pressures whereas
+  Koop 2000 uses the difference in chemical potential. It should be noted that the Koop 2000
+  parameterization is only valid for temperatures up to 240K and a temperature-dependent max
+  pressure. The max valid pressure becomes negative around 237K, so the Koop 2000 parameterizaiton
+  should not be valid beyond 237K.
+
+```@example
+include("HomFreezingPlots.jl")
+```
+![](homogeneous_J.svg)
diff --git a/docs/src/IceNucleationBox.md b/docs/src/IceNucleationBox.md
index ddad8a73bd..206de11c04 100644
--- a/docs/src/IceNucleationBox.md
+++ b/docs/src/IceNucleationBox.md
@@ -4,7 +4,7 @@
 
 This is a 0-dimensional adiabatic parcel model for testing nucleation schemes.
 It is based on [KorolevMazin2003](@cite), as well as the cirrus box model
-  [Karcher2006](@cite), [Tully2022](@cite).
+  [Karcher2006](@cite), [Tully2023](@cite).
 The model solves for saturation in a 0-dimensional
   adiabatic parcel raising with constant velocity.
 
@@ -172,6 +172,6 @@ Between each run the water vapor specific humidity is changed,
 The prescribed vertical velocity is equal to 3.5 cm/s.
 
 ```@example
-include("../../parcel/parcel.jl")
+include("../../parcel/Tully_et_al_2023.jl")
 ```
-![](cirrus_box.svg)
\ No newline at end of file
+![](cirrus_box.svg)
diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl
new file mode 100644
index 0000000000..eb7751a2d6
--- /dev/null
+++ b/parcel/Tully_et_al_2023.jl
@@ -0,0 +1,191 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# boilerplate code to get free parameter values
+include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+
+"""
+    Wrapper for initial condition
+"""
+function get_initial_condition(
+    prs,
+    N_act,
+    p_a,
+    T,
+    q_vap,
+    q_liq,
+    q_ice,
+    N_aerosol,
+)
+    thermo_params = CMP.thermodynamics_params(prs)
+    q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
+    R_a = TD.gas_constant_air(thermo_params, q)
+    R_v = CMP.R_v(prs)
+    e_si = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
+    e = q_vap * p_a * R_v / R_a
+    S_i = e / e_si
+
+    return [S_i, N_act, p_a, T, q_vap, q_ice, N_aerosol]
+end
+
+"""
+    Wrapper for running the simulation following the same framework as in
+    Tully et al 2023 (https://doi.org/10.5194/gmd-16-2957-2023)
+     - the simulation consists of 3 periods mimicking 3 large scale model steps
+     - each period is 30 minutes long
+     - each period is run with user specified constant timestep
+     - the large scale initial conditions between each period are different
+"""
+function run_parcel(FT)
+
+    # Boiler plate code to have access to model parameters and constants
+    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
+    prs = cloud_microphysics_parameters(toml_dict)
+    thermo_params = CMP.thermodynamics_params(prs)
+
+    # Initial conditions for 1st period
+    N_aerosol = FT(2000 * 1e3)
+    N_0 = FT(0)
+    p_0 = FT(20000)
+    T_0 = FT(230)
+    q_vap_0 = FT(0.0003345)
+    q_liq_0 = FT(0)
+    q_ice_0 = FT(0)
+    # Initial conditions for the 2nd period
+    T2 = FT(229.25)
+    q_vap2 = FT(0.00034)
+    # Initial conditions for the 3rd period
+    T3 = FT(228.55)
+    q_vap3 = FT(0.000345)
+
+    # Simulation time
+    t_max = 30 * 60
+
+    # Simulation parameters passed into ODE solver
+    r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles
+    w = FT(3.5 * 1e-2) # updraft speed
+    α_m = FT(0.5) # accomodation coefficient
+    const_dt = 0.1 # model timestep
+    p = (; prs, const_dt, r_nuc, w, α_m)
+
+    # Simulation 1
+    IC1 = get_initial_condition(
+        prs,
+        N_0,
+        p_0,
+        T_0,
+        q_vap_0,
+        q_liq_0,
+        q_ice_0,
+        N_aerosol,
+    )
+    prob1 = ODE.ODEProblem(parcel_model, IC1, (FT(0), t_max), p)
+    sol1 = ODE.solve(
+        prob1,
+        ODE.Euler(),
+        dt = const_dt,
+        reltol = 10 * eps(FT),
+        abstol = 10 * eps(FT),
+    )
+
+    # Simulation 2
+    # (alternatively set T and take q_vap from the previous simulation)
+    #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end])
+    IC2 = get_initial_condition(
+        prs,
+        sol1[2, end],
+        sol1[3, end],
+        sol1[4, end],
+        q_vap2,
+        q_liq_0,
+        sol1[6, end],
+        sol1[7, end],
+    )
+    prob2 = ODE.ODEProblem(
+        parcel_model,
+        IC2,
+        (FT(sol1.t[end]), sol1.t[end] + t_max),
+        p,
+    )
+    sol2 = ODE.solve(
+        prob2,
+        ODE.Euler(),
+        dt = const_dt,
+        reltol = 10 * eps(FT),
+        abstol = 10 * eps(FT),
+    )
+
+    # Simulation 3
+    # (alternatively set T and take q_vap from the previous simulation)
+    #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end])
+    IC3 = get_initial_condition(
+        prs,
+        sol2[2, end],
+        sol2[3, end],
+        sol2[4, end],
+        q_vap3,
+        q_liq_0,
+        sol2[6, end],
+        sol2[7, end],
+    )
+    prob3 = ODE.ODEProblem(
+        parcel_model,
+        IC3,
+        (FT(sol2.t[end]), sol2.t[end] + t_max),
+        p,
+    )
+    sol3 = ODE.solve(
+        prob3,
+        ODE.Euler(),
+        dt = const_dt,
+        reltol = 10 * eps(FT),
+        abstol = 10 * eps(FT),
+    )
+
+    # Plot results
+    fig = MK.Figure(resolution = (800, 600))
+    ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]")
+    ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
+    ax3 = MK.Axis(fig[2, 1], ylabel = "N act [1/dm3]", yscale = log10)
+    ax4 = MK.Axis(fig[2, 2], ylabel = "N areo [1/dm3]")
+    ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Height [m]")
+    ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]", xlabel = "Height [m]")
+
+    MK.ylims!(ax1, 1.0, 1.5)
+    MK.ylims!(ax3, 3, 2e3)
+
+    MK.lines!(ax1, sol1.t * w, sol1[1, :])
+    MK.lines!(ax1, sol2.t * w, sol2[1, :])
+    MK.lines!(ax1, sol3.t * w, sol3[1, :])
+
+    MK.lines!(ax2, sol1.t * w, sol1[4, :])
+    MK.lines!(ax2, sol2.t * w, sol2[4, :])
+    MK.lines!(ax2, sol3.t * w, sol3[4, :])
+
+    MK.lines!(ax3, sol1.t * w, sol1[2, :] * 1e-3)
+    MK.lines!(ax3, sol2.t * w, sol2[2, :] * 1e-3)
+    MK.lines!(ax3, sol3.t * w, sol3[2, :] * 1e-3)
+
+    MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3)
+    MK.lines!(ax4, sol2.t * w, sol2[7, :] * 1e-3)
+    MK.lines!(ax4, sol3.t * w, sol3[7, :] * 1e-3)
+
+    MK.lines!(ax5, sol1.t * w, sol1[5, :] * 1e3)
+    MK.lines!(ax5, sol2.t * w, sol2[5, :] * 1e3)
+    MK.lines!(ax5, sol3.t * w, sol3[5, :] * 1e3)
+
+    MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3)
+    MK.lines!(ax6, sol2.t * w, sol2[6, :] * 1e3)
+    MK.lines!(ax6, sol3.t * w, sol3[6, :] * 1e3)
+
+    MK.save("cirrus_box.svg", fig)
+end
+
+run_parcel(Float64)
diff --git a/parcel/parcel.jl b/parcel/parcel.jl
index c4995a99d6..6c91bbeca6 100644
--- a/parcel/parcel.jl
+++ b/parcel/parcel.jl
@@ -1,9 +1,5 @@
-import OrdinaryDiffEq as ODE
-import CairoMakie as MK
-
 import Thermodynamics as TD
 import CloudMicrophysics as CM
-import CLIMAParameters as CP
 
 const CMT = CM.CommonTypes
 const CMO = CM.Common
@@ -15,7 +11,7 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 """
     ODE problem definition
 """
-function cirrus_box(dY, Y, p, t)
+function parcel_model(dY, Y, p, t)
 
     # Get simulation parameters
     (; prs, const_dt, r_nuc, w, α_m) = p
@@ -99,182 +95,3 @@ function cirrus_box(dY, Y, p, t)
 
     # TODO - add diagnostics output (radius, S, etc)
 end
-
-"""
-    Wrapper for initial condition
-"""
-function get_initial_condition(
-    prs,
-    N_act,
-    p_a,
-    T,
-    q_vap,
-    q_liq,
-    q_ice,
-    N_aerosol,
-)
-    thermo_params = CMP.thermodynamics_params(prs)
-    q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
-    R_a = TD.gas_constant_air(thermo_params, q)
-    R_v = CMP.R_v(prs)
-    e_si = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice())
-    e = q_vap * p_a * R_v / R_a
-    S_i = e / e_si
-
-    return [S_i, N_act, p_a, T, q_vap, q_ice, N_aerosol]
-end
-
-"""
-    Wrapper for running the simulation following the same framework as in
-    Tully et al 2022 (10.5194/egusphere-2022-1057)
-     - the simulation consists of 3 periods mimicking 3 large scale model steps
-     - each period is 30 minutes long
-     - each period is run with user specified constant timestep
-     - the large scale initial conditions between each period are different
-"""
-function run_parcel(FT)
-
-    # Boiler plate code to have access to model parameters and constants
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    prs = cloud_microphysics_parameters(toml_dict)
-    thermo_params = CMP.thermodynamics_params(prs)
-
-    # Initial conditions for 1st period
-    N_aerosol = FT(2000 * 1e3)
-    N_0 = FT(0)
-    p_0 = FT(20000)
-    T_0 = FT(230)
-    q_vap_0 = FT(0.0003345)
-    q_liq_0 = FT(0)
-    q_ice_0 = FT(0)
-    # Initial conditions for the 2nd period
-    T2 = FT(229.25)
-    q_vap2 = FT(0.00034)
-    # Initial conditions for the 3rd period
-    T3 = FT(228.55)
-    q_vap3 = FT(0.000345)
-
-    # Simulation time
-    t_max = 30 * 60
-
-    # Simulation parameters passed into ODE solver
-    r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles
-    w = FT(3.5 * 1e-2) # updraft speed
-    α_m = FT(0.5) # accomodation coefficient
-    const_dt = 0.1 # model timestep
-    p = (; prs, const_dt, r_nuc, w, α_m)
-
-    # Simulation 1
-    IC1 = get_initial_condition(
-        prs,
-        N_0,
-        p_0,
-        T_0,
-        q_vap_0,
-        q_liq_0,
-        q_ice_0,
-        N_aerosol,
-    )
-    prob1 = ODE.ODEProblem(cirrus_box, IC1, (FT(0), t_max), p)
-    sol1 = ODE.solve(
-        prob1,
-        ODE.Euler(),
-        dt = const_dt,
-        reltol = 10 * eps(FT),
-        abstol = 10 * eps(FT),
-    )
-
-    # Simulation 2
-    # (alternatively set T and take q_vap from the previous simulation)
-    #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end])
-    IC2 = get_initial_condition(
-        prs,
-        sol1[2, end],
-        sol1[3, end],
-        sol1[4, end],
-        q_vap2,
-        q_liq_0,
-        sol1[6, end],
-        sol1[7, end],
-    )
-    prob2 = ODE.ODEProblem(
-        cirrus_box,
-        IC2,
-        (FT(sol1.t[end]), sol1.t[end] + t_max),
-        p,
-    )
-    sol2 = ODE.solve(
-        prob2,
-        ODE.Euler(),
-        dt = const_dt,
-        reltol = 10 * eps(FT),
-        abstol = 10 * eps(FT),
-    )
-
-    # Simulation 3
-    # (alternatively set T and take q_vap from the previous simulation)
-    #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end])
-    IC3 = get_initial_condition(
-        prs,
-        sol2[2, end],
-        sol2[3, end],
-        sol2[4, end],
-        q_vap3,
-        q_liq_0,
-        sol2[6, end],
-        sol2[7, end],
-    )
-    prob3 = ODE.ODEProblem(
-        cirrus_box,
-        IC3,
-        (FT(sol2.t[end]), sol2.t[end] + t_max),
-        p,
-    )
-    sol3 = ODE.solve(
-        prob3,
-        ODE.Euler(),
-        dt = const_dt,
-        reltol = 10 * eps(FT),
-        abstol = 10 * eps(FT),
-    )
-
-    # Plot results
-    fig = MK.Figure(resolution = (800, 600))
-    ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]")
-    ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
-    ax3 = MK.Axis(fig[2, 1], ylabel = "N act [1/dm3]", yscale = log10)
-    ax4 = MK.Axis(fig[2, 2], ylabel = "N areo [1/dm3]")
-    ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Height [m]")
-    ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]", xlabel = "Height [m]")
-
-    MK.ylims!(ax1, 1.0, 1.5)
-    MK.ylims!(ax3, 3, 2e3)
-
-    MK.lines!(ax1, sol1.t * w, sol1[1, :])
-    MK.lines!(ax1, sol2.t * w, sol2[1, :])
-    MK.lines!(ax1, sol3.t * w, sol3[1, :])
-
-    MK.lines!(ax2, sol1.t * w, sol1[4, :])
-    MK.lines!(ax2, sol2.t * w, sol2[4, :])
-    MK.lines!(ax2, sol3.t * w, sol3[4, :])
-
-    MK.lines!(ax3, sol1.t * w, sol1[2, :] * 1e-3)
-    MK.lines!(ax3, sol2.t * w, sol2[2, :] * 1e-3)
-    MK.lines!(ax3, sol3.t * w, sol3[2, :] * 1e-3)
-
-    MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3)
-    MK.lines!(ax4, sol2.t * w, sol2[7, :] * 1e-3)
-    MK.lines!(ax4, sol3.t * w, sol3[7, :] * 1e-3)
-
-    MK.lines!(ax5, sol1.t * w, sol1[5, :] * 1e3)
-    MK.lines!(ax5, sol2.t * w, sol2[5, :] * 1e3)
-    MK.lines!(ax5, sol3.t * w, sol3[5, :] * 1e3)
-
-    MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3)
-    MK.lines!(ax6, sol2.t * w, sol2[6, :] * 1e3)
-    MK.lines!(ax6, sol3.t * w, sol3[6, :] * 1e3)
-
-    MK.save("cirrus_box.svg", fig)
-end
-
-run_parcel(Float64)
diff --git a/src/Common.jl b/src/Common.jl
index 01f3523be5..9f7f041d77 100644
--- a/src/Common.jl
+++ b/src/Common.jl
@@ -129,7 +129,7 @@ Returns the saturation vapor pressure above a sulphuric acid solution droplet in
 """
 function H2SO4_soln_saturation_vapor_pressure(x::FT, T::FT) where {FT <: Real}
 
-    @assert T < FT(235)
+    @assert T < FT(240)
     @assert T > FT(185)
 
     w_h = 1.4408 * x
diff --git a/test/common_functions_tests.jl b/test/common_functions_tests.jl
index c78105e6a6..24af5dbfaa 100644
--- a/test/common_functions_tests.jl
+++ b/test/common_functions_tests.jl
@@ -58,12 +58,12 @@ function test_H2SO4_soln_saturation_vapor_pressure(FT)
 
         T_warm = FT(225.0)
         T_cold = FT(200.0)
-        T_too_warm = FT(240)
+        T_too_warm = FT(245)
         T_too_cold = FT(180)
         x_sulph = FT(0.1)
 
         # If T out of range
-        TT.@test_throws AssertionError("T < FT(235)") CO.H2SO4_soln_saturation_vapor_pressure(
+        TT.@test_throws AssertionError("T < FT(240)") CO.H2SO4_soln_saturation_vapor_pressure(
             x_sulph,
             T_too_warm,
         )
diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl
index 9c595253e2..e7567a36d6 100644
--- a/test/heterogeneous_ice_nucleation_tests.jl
+++ b/test/heterogeneous_ice_nucleation_tests.jl
@@ -6,7 +6,7 @@ import CLIMAParameters as CP
 
 const CMT = CM.CommonTypes
 const CO = CM.Common
-const CMI = CM.HetIceNucleation
+const CMI_het = CM.HetIceNucleation
 const ArizonaTestDust = CMT.ArizonaTestDustType()
 const DesertDust = CMT.DesertDustType()
 
@@ -14,7 +14,8 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 
 @info "Heterogeneous Ice Nucleation Tests"
 
-function test_dust_activation(FT)
+function test_heterogeneous_ice_nucleation(FT)
+
     toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
     prs = cloud_microphysics_parameters(toml_dict)
 
@@ -30,7 +31,7 @@ function test_dust_activation(FT)
         # No activation below critical supersaturation
         for dust in [ArizonaTestDust, DesertDust]
             for T in [T_warm, T_cold]
-                TT.@test CMI.dust_activated_number_fraction(
+                TT.@test CMI_het.dust_activated_number_fraction(
                     prs,
                     Si_low,
                     T,
@@ -41,23 +42,23 @@ function test_dust_activation(FT)
 
         # Activate more in cold temperatures and higher supersaturations
         for dust in [ArizonaTestDust, DesertDust]
-            TT.@test CMI.dust_activated_number_fraction(
+            TT.@test CMI_het.dust_activated_number_fraction(
                 prs,
                 Si_hgh,
                 T_warm,
                 dust,
-            ) > CMI.dust_activated_number_fraction(
+            ) > CMI_het.dust_activated_number_fraction(
                 prs,
                 Si_med,
                 T_warm,
                 dust,
             )
-            TT.@test CMI.dust_activated_number_fraction(
+            TT.@test CMI_het.dust_activated_number_fraction(
                 prs,
                 Si_med,
                 T_cold,
                 dust,
-            ) > CMI.dust_activated_number_fraction(
+            ) > CMI_het.dust_activated_number_fraction(
                 prs,
                 Si_med,
                 T_warm,
@@ -67,7 +68,7 @@ function test_dust_activation(FT)
 
         for dust in [ArizonaTestDust, DesertDust]
             for T in [T_warm, T_cold]
-                TT.@test CMI.dust_activated_number_fraction(
+                TT.@test CMI_het.dust_activated_number_fraction(
                     prs,
                     Si_too_hgh,
                     T,
@@ -76,11 +77,6 @@ function test_dust_activation(FT)
             end
         end
     end
-end
-
-function test_ABIFM_J(FT)
-    toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
-    prs = cloud_microphysics_parameters(toml_dict)
 
     TT.@testset "ABIFM J" begin
 
@@ -91,25 +87,15 @@ function test_ABIFM_J(FT)
         # higher nucleation rate at colder temperatures
         for dust in
             [CMT.IlliteType(), CMT.KaoliniteType(), CMT.DesertDustType()]
-            TT.@test CMI.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_cold)) >
-                     CMI.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_warm))
+            TT.@test CMI_het.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_cold)) >
+                     CMI_het.ABIFM_J(dust, CO.Delta_a_w(prs, x_sulph, T_warm))
         end
     end
 end
 
-
-
-println("Testing Float64")
-test_dust_activation(Float64)
-
-
-println("Testing Float32")
-test_dust_activation(Float32)
-
-
 println("Testing Float64")
-test_ABIFM_J(Float64)
+test_heterogeneous_ice_nucleation(Float64)
 
 
 println("Testing Float32")
-test_ABIFM_J(Float32)
+test_heterogeneous_ice_nucleation(Float32)
diff --git a/test/homogeneous_ice_nucleation_tests.jl b/test/homogeneous_ice_nucleation_tests.jl
index 5a06fea86f..676b7fd6f0 100644
--- a/test/homogeneous_ice_nucleation_tests.jl
+++ b/test/homogeneous_ice_nucleation_tests.jl
@@ -5,7 +5,7 @@ import CloudMicrophysics as CM
 import CLIMAParameters as CP
 
 const CO = CM.Common
-const CMI = CM.HomIceNucleation
+const CMI_hom = CM.HomIceNucleation
 
 include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 
@@ -24,14 +24,14 @@ function test_homogeneous_J(FT)
         Δa_w_too_large = FT(0.35)
 
         # higher nucleation rate at colder temperatures
-        TT.@test CMI.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_cold)) >
-                 CMI.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_warm))
+        TT.@test CMI_hom.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_cold)) >
+                 CMI_hom.homogeneous_J(CO.Delta_a_w(prs, x_sulph, T_warm))
 
         # If Δa_w out of range
-        TT.@test_throws AssertionError("Δa_w > 0.26") CMI.homogeneous_J(
+        TT.@test_throws AssertionError("Δa_w > 0.26") CMI_hom.homogeneous_J(
             Δa_w_too_small,
         )
-        TT.@test_throws AssertionError("Δa_w < 0.34") CMI.homogeneous_J(
+        TT.@test_throws AssertionError("Δa_w < 0.34") CMI_hom.homogeneous_J(
             Δa_w_too_large,
         )
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index ebe77e97f6..e3885b252e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 include("performance_tests.jl")
 include("aerosol_activation_tests.jl")
 include("heterogeneous_ice_nucleation_tests.jl")
+include("homogeneous_ice_nucleation_tests.jl")
 include("microphysics_tests.jl")
 include("gpu_tests.jl")
 include("common_functions_tests.jl")