diff --git a/CHANGELOG.md b/CHANGELOG.md
index dacfb37..1eb1bd2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,17 +2,25 @@
 
 All notable changes to this project will be documented in this file.
 
-## [Unreleased]
-
+## [v0.10.0] 2024-11-19
+
+- **ADDED** Top-level function exports to avoid submodule imports
+  - `pymsis.calculate()` is the primrary entrypoint to running the MSIS
+    model and calculating the atmosphere at the requested data points.
+  - `pymsis.msis.run()` was not very descriptive and caused issues with IDL
+    bridging into Python and wanting that name reserved. To avoid this, the
+    function has been renamed `calculate` and is available as `pymsis.calculate()`
+    now. The `pymsis.msis.run()` is still available as an alias for now, but
+    may be removed in the future.
+- **ADDED** Variable enumeration for easier indexing into output arrays.
+  - This can be used as `pymsis.Variable.O2` for getting the `O2` species index.
+    For example, `output_array[..., pymsis.Variable.O2]`.
 - **ADDED** Python 3.13 and 3.13t support
 - **ADDED** Multithreaded support.
   - The underlying MSIS libraries are not threadsafe due
     to the use of many global/save variables. There is a lock around the
     extension modules so that only one thread will be calling the routines
     at a time, so the Python library is safe to use in a multi-threaded context.
-- **ADDED** Variable enumeration for easier indexing into output arrays.
-  - This can be used as `msis.Variable.O2` for getting the `O2` species index.
-    For example, `output_array[..., msis.Variable.O2]`.
 - **MAINTENANCE** Default `-O1` optimization level for all builds.
   - Previously, this
     was only done on Windows machines. Users can change this by updating
@@ -30,7 +38,7 @@ All notable changes to this project will be documented in this file.
 - **DEPRECATED** Calling `msis00f.pytselec()` and `msis00f.pygtd7d` functions.
   - Use `msis00f.pyinitswitch` and `msis00f.pymsiscalc` instead.
   - This helps with standardization across the extension modules. These extension
-    should rarely be used by external people and `msis.run()` is a better entry
+    should rarely be used by external people and `pymsis.calculate()` is a better entry
     to using the package.
 
 ## [v0.9.0] - 2024-04-03
@@ -57,7 +65,7 @@ All notable changes to this project will be documented in this file.
 ## [v0.6.0] - 2022-11-14
 
 - **ADDED** Automatic download of F10.7 and ap data for users.
-  - This means that F10.7 and ap are optional inputs to the `msis.run()`
+  - This means that F10.7 and ap are optional inputs to the `pymsis.calculate()`
     function during historical periods and the routines will automatically
     fetch the proper input data.
 
@@ -67,7 +75,7 @@ All notable changes to this project will be documented in this file.
 
 - **ADDED** MSIS2.1, a new version of MSIS.
   - This is the first version that contains NO.
-  - This is the new default used in `msis.run()`.
+  - This is the new default used in `pymsis.calculate()`.
 - **MAINTENANCE** Added more wheels to the release and CI systems for testing.
 
 ## [v0.4.0] - 2022-02-26
diff --git a/README.md b/README.md
index 35023d0..3247dec 100644
--- a/README.md
+++ b/README.md
@@ -28,11 +28,11 @@ there is a web viewer that uses pymsis: <https://swx-trec.com/msis>
 
 ```python
 import numpy as np
-from pymsis import msis
+import pymsis
 
 dates = np.arange(np.datetime64("2003-10-28T00:00"), np.datetime64("2003-11-04T00:00"), np.timedelta64(30, "m"))
 # geomagnetic_activity=-1 is a storm-time run
-data = msis.run(dates, 0, 0, 400, geomagnetic_activity=-1)
+data = pymsis.calculate(dates, 0, 0, 400, geomagnetic_activity=-1)
 
 # Plot the data
 import matplotlib.pyplot as plt
@@ -56,7 +56,7 @@ is developed by the Naval Research Laboratory.
 Note that the MSIS2 code is not available for commercial use without
 contacting NRL. See the [MSIS2 license file](https://github.com/SWxTREC/pymsis/blob/main/MSIS2_LICENSE)) for explicit
 details. We do not repackage the MSIS source code in this
-repository for that reason. However, we do provide utilities to easily
+repository for that reason. However, utility functions are provided to easily
 download and extract the original source code. By using that code you
 agree to their terms and conditions.
 
diff --git a/docs/source/reference/index.rst b/docs/source/reference/index.rst
index f1e18da..07b0571 100644
--- a/docs/source/reference/index.rst
+++ b/docs/source/reference/index.rst
@@ -5,30 +5,47 @@ API reference
 .. currentmodule:: pymsis
 
 This page gives an overview of the routines within the pymsis package.
-To run the code, use the :mod:`pymsis.msis` module.
+To calculate atmospheric constituents at grid points, use the :func:`~.calculate` function.
+This is typically the only function you will need to use. You can access the output variables
+using the :class:`~.Variable` enumeration for easier indexing into the output data array.
+For example, as ``output_array[..., Variable.MASS_DENSITY]``.
 
 .. autosummary::
     :toctree: generated/
+    :nosignatures:
 
-    msis.run
-    msis.create_options
-    msis.create_input
+    calculate
+    Variable
+
+msis module
+-----------
 
-The output of the model is stored in basic numpy arrays with the final
-dimension being the variable/species. To get the output in a more human-readable
-format, use the :class:`~.Variable` enumeration that provides
-easier indexing into the output arrays.
+For more control and help creating properly formatted inputs, one can
+use the :mod:`pymsis.msis` module. This module provides functions to
+create input data, create the options list, and run the model.
 
 .. autosummary::
     :toctree: generated/
     :nosignatures:
 
-    msis.Variable
+    msis.create_input
+    msis.create_options
+    msis.calculate
 
-To get input data for historical events, use the :mod:`pymsis.utils` module.
+utils module
+------------
+
+The MSIS model requires solar forcing data (F10.7 and ap) to calculate the state of the atmosphere.
+When running over past time periods, this data is automatically downloaded and used if not specified by the user.
+For more control over the F10.7 and Ap values the model is using, one can input the desired values
+into the :func:`~.calculate` call ``calculate(..., f107s=..., f107as=..., aps=...)``.
+To get the values used automatically by the model, the :mod:`pymsis.utils` module has several
+helper routines to download new files (will retrieve a file with data up to the present) and
+get the proper F10.7 and ap values to use for a specific date and time.
 
 .. autosummary::
     :toctree: generated/
+    :nosignatures:
 
     utils.download_f107_ap
-    utils.get_f107_ap
\ No newline at end of file
+    utils.get_f107_ap
diff --git a/examples/plot_altitude_profiles.py b/examples/plot_altitude_profiles.py
index 479e84a..a3edcc6 100644
--- a/examples/plot_altitude_profiles.py
+++ b/examples/plot_altitude_profiles.py
@@ -12,7 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lon = 0
@@ -24,9 +24,9 @@
 aps = [[ap] * 7]
 
 date = np.datetime64("2003-01-01T00:00")
-output_midnight = msis.run(date, lon, lat, alts, f107, f107a, aps)
+output_midnight = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps)
 date = np.datetime64("2003-01-01T12:00")
-output_noon = msis.run(date, lon, lat, alts, f107, f107a, aps)
+output_noon = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps)
 
 #  output is now of the shape (1, 1, 1, 1000, 11)
 # Get rid of the single dimensions
@@ -34,7 +34,7 @@
 output_noon = np.squeeze(output_noon)
 
 _, ax = plt.subplots()
-for variable in msis.Variable:
+for variable in pymsis.Variable:
     if variable.name in ("Total mass density", "Temperature"):
         # Ignore non-number densities
         continue
diff --git a/examples/plot_annual_variation.py b/examples/plot_annual_variation.py
index ea7e3ea..a4c1a92 100644
--- a/examples/plot_annual_variation.py
+++ b/examples/plot_annual_variation.py
@@ -10,7 +10,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lon = 0
@@ -27,7 +27,7 @@
 f107as = [f107a] * ndates
 aps = [[ap] * 7] * ndates
 
-output = msis.run(dates, lon, lat, alt, f107s, f107as, aps)
+output = pymsis.calculate(dates, lon, lat, alt, f107s, f107as, aps)
 #  output is now of the shape (ndates, 1, 1, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
@@ -36,7 +36,7 @@
 variation = 100 * (output / output.mean(axis=0) - 1)
 
 _, ax = plt.subplots()
-for variable in msis.Variable:
+for variable in pymsis.Variable:
     if variable.name == "NO":
         # There is currently no NO data
         continue
diff --git a/examples/plot_diurnal_variation.py b/examples/plot_diurnal_variation.py
index eee56cb..e002bd6 100644
--- a/examples/plot_diurnal_variation.py
+++ b/examples/plot_diurnal_variation.py
@@ -11,7 +11,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lon = 0
@@ -28,7 +28,7 @@
 f107as = [f107a] * ndates
 aps = [[ap] * 7] * ndates
 
-output = msis.run(dates, lon, lat, alt, f107s, f107as, aps)
+output = pymsis.calculate(dates, lon, lat, alt, f107s, f107as, aps)
 #  output is now of the shape (ndates, 1, 1, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
@@ -37,7 +37,7 @@
 variation = 100 * (output / output.mean(axis=0) - 1)
 
 _, ax = plt.subplots()
-for variable in msis.Variable:
+for variable in pymsis.Variable:
     ax.plot(dates, variation[:, variable], label=variable.name)
 
 ax.legend(
diff --git a/examples/plot_surface.py b/examples/plot_surface.py
index d42fb67..8f62495 100644
--- a/examples/plot_surface.py
+++ b/examples/plot_surface.py
@@ -10,7 +10,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lons = range(-180, 185, 5)
@@ -23,14 +23,14 @@
 date = np.datetime64("2003-01-01T12:00")
 aps = [[ap] * 7]
 
-output = msis.run(date, lons, lats, alt, f107, f107a, aps)
+output = pymsis.calculate(date, lons, lats, alt, f107, f107a, aps)
 #  output is now of the shape (1, nlons, nlats, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
 
 _, ax = plt.subplots()
 xx, yy = np.meshgrid(lons, lats)
-mesh = ax.pcolormesh(xx, yy, output[:, :, msis.Variable.O2].T, shading="auto")
+mesh = ax.pcolormesh(xx, yy, output[:, :, pymsis.Variable.O2].T, shading="auto")
 plt.colorbar(mesh, label="Number density (/m$^3$)")
 
 ax.set_title(f"O$_2$ number density at {alt} km")
diff --git a/examples/plot_surface_animation.py b/examples/plot_surface_animation.py
index 913fdfa..d57a050 100644
--- a/examples/plot_surface_animation.py
+++ b/examples/plot_surface_animation.py
@@ -13,7 +13,7 @@
 import numpy as np
 from matplotlib.animation import FuncAnimation
 
-from pymsis import msis
+import pymsis
 
 
 lons = range(-180, 185, 5)
@@ -30,7 +30,7 @@
 f107as = [f107a] * ndates
 aps = [[ap] * 7] * ndates
 
-output = msis.run(dates, lons, lats, alt, f107s, f107as, aps)
+output = pymsis.calculate(dates, lons, lats, alt, f107s, f107as, aps)
 #  output is now of the shape (ndates, nlons, nlats, 1, 11)
 # Get rid of the single dimensions
 output = np.squeeze(output)
@@ -42,14 +42,14 @@
 vmin, vmax = 1e13, 8e13
 norm = matplotlib.colors.Normalize(vmin, vmax)
 mesh = ax_mesh.pcolormesh(
-    xx, yy, output[0, :, :, msis.Variable.N].T, shading="auto", norm=norm
+    xx, yy, output[0, :, :, pymsis.Variable.N].T, shading="auto", norm=norm
 )
 plt.colorbar(
     mesh, label=f"N number density at {alt} km (/m$^3$)", orientation="horizontal"
 )
 
 time_data = output[:, len(lons) // 2, len(lats) // 2, :]
-ax_time.plot(dates, time_data[:, msis.Variable.N], c="k")
+ax_time.plot(dates, time_data[:, pymsis.Variable.N], c="k")
 ax_time.set_xlim(dates[0], dates[-1])
 ax_time.set_ylim(vmin, vmax)
 ax_time.xaxis.set_major_locator(mdates.HourLocator(interval=3))
@@ -74,7 +74,7 @@ def update_surface(t):
     date_string = dates[t].astype("O").strftime("%H:%M")
     title.set_text(f"{date_string}")
     time_line.set_xdata([dates[t]])
-    mesh.set_array(output[t, :, :, msis.Variable.N].T)
+    mesh.set_array(output[t, :, :, pymsis.Variable.N].T)
     sun.set_data([sun_loc[t]], [0])
 
 
diff --git a/examples/plot_version_diff_altitude.py b/examples/plot_version_diff_altitude.py
index b2141ce..fc6d9b2 100644
--- a/examples/plot_version_diff_altitude.py
+++ b/examples/plot_version_diff_altitude.py
@@ -12,7 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lon = 0
@@ -24,13 +24,13 @@
 aps = [[ap] * 7]
 
 date = np.datetime64("2003-01-01T00:00")
-output_midnight2 = msis.run(date, lon, lat, alts, f107, f107a, aps)
-output_midnight0 = msis.run(date, lon, lat, alts, f107, f107a, aps, version=0)
+output_midnight2 = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps)
+output_midnight0 = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps, version=0)
 diff_midnight = (output_midnight2 - output_midnight0) / output_midnight0 * 100
 
 date = np.datetime64("2003-01-01T12:00")
-output_noon2 = msis.run(date, lon, lat, alts, f107, f107a, aps)
-output_noon0 = msis.run(date, lon, lat, alts, f107, f107a, aps, version=0)
+output_noon2 = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps)
+output_noon0 = pymsis.calculate(date, lon, lat, alts, f107, f107a, aps, version=0)
 diff_noon = (output_noon2 - output_noon0) / output_noon0 * 100
 
 
@@ -40,7 +40,7 @@
 diff_noon = np.squeeze(diff_noon)
 
 _, ax = plt.subplots()
-for variable in msis.Variable:
+for variable in pymsis.Variable:
     if variable.name in ("NO", "Total mass density", "Temperature"):
         # There is currently no NO data for earlier versions,
         # also ignore non-number densities
diff --git a/examples/plot_version_diff_surface.py b/examples/plot_version_diff_surface.py
index dc2aad3..651e108 100644
--- a/examples/plot_version_diff_surface.py
+++ b/examples/plot_version_diff_surface.py
@@ -12,7 +12,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pymsis import msis
+import pymsis
 
 
 lons = range(-180, 185, 5)
@@ -25,8 +25,8 @@
 date = np.datetime64("2003-01-01T12:00")
 aps = [[ap] * 7]
 
-output2 = msis.run(date, lons, lats, alt, f107, f107a, aps)
-output0 = msis.run(date, lons, lats, alt, f107, f107a, aps, version=0)
+output2 = pymsis.calculate(date, lons, lats, alt, f107, f107a, aps)
+output0 = pymsis.calculate(date, lons, lats, alt, f107, f107a, aps, version=0)
 diff = (output2 - output0) / output0 * 100
 #  diff is now of the shape (1, nlons, nlats, 1, 11)
 # Get rid of the single dimensions
@@ -36,7 +36,7 @@
 xx, yy = np.meshgrid(lons, lats)
 norm = mpl.colors.Normalize(-50, 50)
 cmap = mpl.colormaps["RdBu_r"]
-for i, variable in enumerate(msis.Variable):
+for i, variable in enumerate(pymsis.Variable):
     if i > 8:
         break
     ax = axarr.flatten()[i]
diff --git a/pymsis/__init__.py b/pymsis/__init__.py
index c99d599..05d45fa 100644
--- a/pymsis/__init__.py
+++ b/pymsis/__init__.py
@@ -2,5 +2,9 @@
 
 import importlib.metadata
 
+from pymsis.msis import Variable, calculate
+
 
 __version__ = importlib.metadata.version("pymsis")
+
+__all__ = ["__version__", "Variable", "calculate"]
diff --git a/pymsis/msis.py b/pymsis/msis.py
index 6ac0679..a75205a 100644
--- a/pymsis/msis.py
+++ b/pymsis/msis.py
@@ -23,8 +23,8 @@
 
 
 class Variable(IntEnum):
-    """
-    Enumeration of output data indices for the pymsis run calls.
+    r"""
+    Enumeration of variable data indices for the output from ``calculate()``.
 
     This can be used to access the data from the output arrays instead of having
     to remember the order of the output. For example,
@@ -33,25 +33,25 @@ class Variable(IntEnum):
     Attributes
     ----------
     MASS_DENSITY
-        Index of total mass density (kg/m3).
+        Index of total mass density (kg/m\ :sup:`3`).
     N2
-        Index of N2 number density (m-3).
+        Index of N2 number density (m\ :sup:`-3`).
     O2
-        Index of O2 number density (m-3).
+        Index of O2 number density (m\ :sup:`-3`).
     O
-        Index of O number density (m-3).
+        Index of O number density (m\ :sup:`-3`).
     HE
-        Index of He number density (m-3).
+        Index of He number density (m\ :sup:`-3`).
     H
-        Index of H number density (m-3).
+        Index of H number density (m\ :sup:`-3`).
     AR
-        Index of Ar number density (m-3).
+        Index of Ar number density (m\ :sup:`-3`).
     N
-        Index of N number density (m-3).
+        Index of N number density (m\ :sup:`-3`).
     ANOMALOUS_O
-        Index of anomalous oxygen number density (m-3).
+        Index of anomalous oxygen number density (m\ :sup:`-3`).
     NO
-        Index of NO number density (m-3).
+        Index of NO number density (m\ :sup:`-3`).
     TEMPERATURE
         Index of temperature (K).
 
@@ -70,7 +70,7 @@ class Variable(IntEnum):
     TEMPERATURE = 10
 
 
-def run(
+def calculate(
     dates: npt.ArrayLike,
     lons: npt.ArrayLike,
     lats: npt.ArrayLike,
@@ -83,7 +83,7 @@ def run(
     version: float | str = 2.1,
     **kwargs: dict,
 ) -> npt.NDArray:
-    """
+    r"""
     Call MSIS to calculate the atmosphere at the provided input points.
 
     **Satellite Fly-Through Mode:**
@@ -136,23 +136,23 @@ def run(
         MSIS version number, one of (0, 2.0, 2.1).
     **kwargs : dict
         Single options for the switches can be defined through keyword arguments.
-        For example, run(..., geomagnetic_activity=-1) will set the geomagnetic
+        For example, calculate(..., geomagnetic_activity=-1) will set the geomagnetic
         activity switch to -1 (storm-time ap mode).
 
     Returns
     -------
     ndarray (ndates, nlons, nlats, nalts, 11) or (ndates, 11)
         | The data calculated at each grid point:
-        | [Total mass density (kg/m3),
-        | N2 # density (m-3),
-        | O2 # density (m-3),
-        | O # density (m-3),
-        | He # density (m-3),
-        | H # density (m-3),
-        | Ar # density (m-3),
-        | N # density (m-3),
-        | Anomalous oxygen # density (m-3),
-        | NO # density (m-3),
+        | [Total mass density (kg/m\ :sup:`3`),
+        | N2 # density (m\ :sup:`-3`),
+        | O2 # density (m\ :sup:`-3`),
+        | O # density (m\ :sup:`-3`),
+        | He # density (m\ :sup:`-3`),
+        | H # density (m\ :sup:`-3`),
+        | Ar # density (m\ :sup:`-3`),
+        | N # density (m\ :sup:`-3`),
+        | Anomalous oxygen # density (m\ :sup:`-3`),
+        | NO # density (m\ :sup:`-3`),
         | Temperature (K)]
 
     Other Parameters
@@ -248,6 +248,10 @@ def run(
     return output.reshape(*input_shape, 11)
 
 
+# For backwards compatibility export the old name here
+run = calculate
+
+
 def create_options(
     f107: float = 1,
     time_independent: float = 1,
diff --git a/pymsis/utils.py b/pymsis/utils.py
index 72f060f..f2f2797 100644
--- a/pymsis/utils.py
+++ b/pymsis/utils.py
@@ -22,13 +22,13 @@ def download_f107_ap() -> None:
     Download the latest ap and F10.7 values.
 
     The file is downloaded into the installed package location, keeping
-    the same filename as the data source: `SW-All.csv`.
+    the same filename as the data source: ``SW-All.csv``.
     This routine can be called to update the data as well if you would like to
     use newer data since the last time you downloaded the file.
 
     Notes
     -----
-    The data provider we are using is CelesTrak, which gets data from other
+    The default data provider is CelesTrak, which gets data from other
     sources and interpolates or predicts missing values to make data easier to
     work with. A warning is emitted when the interpolation or prediction
     occurs, but it is up to the user to verify the Ap and F10.7 data that is used
@@ -192,7 +192,7 @@ def get_f107_ap(dates: npt.ArrayLike) -> tuple[npt.NDArray, npt.NDArray, npt.NDA
         f107a : np.ndarray
             F10.7 running 81-day average centered on the given date(s)
         ap : np.ndarray
-            | Ap for the given date(s), (1-6 only used if `geomagnetic_activity=-1`)
+            | Ap for the given date(s), (1-6 only used if ``geomagnetic_activity=-1``)
             | (0) Daily Ap
             | (1) 3 hr ap index for current time
             | (2) 3 hr ap index for 3 hrs before current time
diff --git a/tests/test_msis.py b/tests/test_msis.py
index 3b9121b..ce6dec4 100644
--- a/tests/test_msis.py
+++ b/tests/test_msis.py
@@ -5,6 +5,7 @@
 import pytest
 from numpy.testing import assert_allclose, assert_array_equal
 
+import pymsis
 from pymsis import msis, msis00f, msis20f, msis21f
 
 
@@ -221,41 +222,47 @@ def test_create_input_lon_wrapping(input_data, expected_input, version):
 
     # Test that -90 and 270 produce the same output
     assert_array_equal(
-        msis.run(date, lons, [lat] * 5, alt, f107, f107a, ap, version=version),
-        msis.run(date, lons + 360, [lat] * 5, alt, f107, f107a, ap, version=version),
+        pymsis.calculate(date, lons, [lat] * 5, alt, f107, f107a, ap, version=version),
+        pymsis.calculate(
+            date, lons + 360, [lat] * 5, alt, f107, f107a, ap, version=version
+        ),
     )
 
 
-def test_run_options(input_data, expected_output):
+def test_calculate_options(input_data, expected_output):
     # Default options is all 1's, so make sure they are equivalent
     assert_allclose(
-        np.squeeze(msis.run(*input_data, options=None)), expected_output, rtol=1e-5
+        np.squeeze(pymsis.calculate(*input_data, options=None)),
+        expected_output,
+        rtol=1e-5,
     )
     assert_allclose(
-        np.squeeze(msis.run(*input_data, options=[1] * 25)), expected_output, rtol=1e-5
+        np.squeeze(pymsis.calculate(*input_data, options=[1] * 25)),
+        expected_output,
+        rtol=1e-5,
     )
 
     with pytest.raises(ValueError, match="options needs to be a list"):
-        msis.run(*input_data, options=[1] * 22)
+        pymsis.calculate(*input_data, options=[1] * 22)
 
 
-def test_run_single_point(input_data, expected_output):
-    output = msis.run(*input_data)
+def test_calculate_single_point(input_data, expected_output):
+    output = pymsis.calculate(*input_data)
     assert output.shape == (1, 11)
     assert_allclose(np.squeeze(output), expected_output, rtol=1e-5)
 
 
-def test_run_gridded_multi_point(input_data, expected_output):
+def test_calculate_gridded_multi_point(input_data, expected_output):
     date, lon, lat, alt, f107, f107a, ap = input_data
     # 5 x 5 surface
     input_data = (date, [lon] * 5, [lat] * 5, alt, f107, f107a, ap)
-    output = msis.run(*input_data)
+    output = pymsis.calculate(*input_data)
     assert output.shape == (1, 5, 5, 1, 11)
     expected = np.tile(expected_output, (5, 5, 1))
     assert_allclose(np.squeeze(output), expected, rtol=1e-5)
 
 
-def test_run_multi_point(input_data, expected_output):
+def test_calculate_multi_point(input_data, expected_output):
     # test multi-point run, like a satellite fly-through
     # and make sure we don't grid the input data
     # 5 input points
@@ -269,77 +276,77 @@ def test_run_multi_point(input_data, expected_output):
         [f107a] * 5,
         ap * 5,
     )
-    output = msis.run(*input_data)
+    output = pymsis.calculate(*input_data)
     assert output.shape == (5, 11)
     expected = np.tile(expected_output, (5, 1))
     assert_allclose(np.squeeze(output), expected, rtol=1e-5)
 
 
-def test_run_wrapped_lon(input_data, expected_output):
+def test_calculate_wrapped_lon(input_data, expected_output):
     date, _, lat, alt, f107, f107a, ap = input_data
 
     input_data = (date, -180, lat, alt, f107, f107a, ap)
-    output1 = msis.run(*input_data)
+    output1 = pymsis.calculate(*input_data)
     input_data = (date, 180, lat, alt, f107, f107a, ap)
-    output2 = msis.run(*input_data)
+    output2 = pymsis.calculate(*input_data)
     assert_allclose(output1, output2, rtol=1e-5)
 
     input_data = (date, 0, lat, alt, f107, f107a, ap)
-    output1 = msis.run(*input_data)
+    output1 = pymsis.calculate(*input_data)
     input_data = (date, 360, lat, alt, f107, f107a, ap)
-    output2 = msis.run(*input_data)
+    output2 = pymsis.calculate(*input_data)
     assert_allclose(output1, output2, rtol=1e-5)
 
 
-def test_run_poles(input_data):
+def test_calculate_poles(input_data):
     # Test that moving in longitude around a pole
     # returns the same values
     # North pole
     date, _, _, alt, f107, f107a, ap = input_data
     input_data = (date, 0, 90, alt, f107, f107a, ap)
-    output1 = msis.run(*input_data)
+    output1 = pymsis.calculate(*input_data)
     input_data = (date, 20, 90, alt, f107, f107a, ap)
-    output2 = msis.run(*input_data)
+    output2 = pymsis.calculate(*input_data)
     assert_allclose(output1, output2, rtol=1e-5)
 
     # South pole
     input_data = (date, 0, -90, alt, f107, f107a, ap)
-    output1 = msis.run(*input_data)
+    output1 = pymsis.calculate(*input_data)
     input_data = (date, 20, -90, alt, f107, f107a, ap)
-    output2 = msis.run(*input_data)
+    output2 = pymsis.calculate(*input_data)
     assert_allclose(output1, output2, rtol=1e-5)
 
 
-def test_run_versions(input_data):
+def test_calculate_versions(input_data):
     # Make sure we accept these version numbers and don't throw an error
     for ver in [0, 2, 2.0, 2.1, "00", "0", "2.0", "2.1", "2"]:
-        msis.run(*input_data, version=ver)
+        pymsis.calculate(*input_data, version=ver)
     # Test for something outside of that list for an error to be thrown
     with pytest.raises(ValueError, match="The MSIS version selected"):
-        msis.run(*input_data, version=1)
+        pymsis.calculate(*input_data, version=1)
 
 
-def test_run_version00(input_data, expected_output00):
-    output = msis.run(*input_data, version=0)
+def test_calculate_version00(input_data, expected_output00):
+    output = pymsis.calculate(*input_data, version=0)
     assert output.shape == (1, 11)
     assert_allclose(np.squeeze(output), expected_output00, rtol=1e-5)
 
 
-def test_run_multi_point00(input_data, expected_output00):
+def test_calculate_multi_point00(input_data, expected_output00):
     date, lon, lat, alt, f107, f107a, ap = input_data
     # 5 x 5 surface
     input_data = (date, [lon] * 5, [lat] * 5, alt, f107, f107a, ap)
-    output = msis.run(*input_data, version=0)
+    output = pymsis.calculate(*input_data, version=0)
     assert output.shape == (1, 5, 5, 1, 11)
     expected = np.tile(expected_output00, (5, 5, 1))
     assert_allclose(np.squeeze(output), expected, rtol=1e-5)
 
 
-def test_run_version00_low_altitude(input_data):
+def test_calculate_version00_low_altitude(input_data):
     # There is no O, H, N below 72.5 km, should be NaN
     date, lon, lat, _, f107, f107a, ap = input_data
     input_data = (date, lon, lat, 71, f107, f107a, ap)
-    output = msis.run(*input_data, version=0)
+    output = pymsis.calculate(*input_data, version=0)
     assert np.all(np.isnan(np.squeeze(output)[[3, 5, 7]]))
 
 
@@ -369,25 +376,25 @@ def test_create_input_auto_f107(input_auto_f107_ap):
     assert_allclose(data, test_data)
 
 
-def test_run_auto_f107(input_auto_f107_ap):
+def test_calculate_auto_f107(input_auto_f107_ap):
     # Dropping any of the f107/f107a/ap will go and get those
     # values automatically as needed
     date, lon, lat, alt, f107, f107a, ap = input_auto_f107_ap
 
     # fully specified run
-    expected = msis.run(*input_auto_f107_ap)
+    expected = pymsis.calculate(*input_auto_f107_ap)
 
     # auto
-    assert_allclose(msis.run(date, lon, lat, alt), expected)
+    assert_allclose(pymsis.calculate(date, lon, lat, alt), expected)
 
     # f107
-    assert_allclose(msis.run(date, lon, lat, alt, f107s=f107), expected)
+    assert_allclose(pymsis.calculate(date, lon, lat, alt, f107s=f107), expected)
 
     # f107a
-    assert_allclose(msis.run(date, lon, lat, alt, f107as=f107a), expected)
+    assert_allclose(pymsis.calculate(date, lon, lat, alt, f107as=f107a), expected)
 
     # ap
-    assert_allclose(msis.run(date, lon, lat, alt, aps=ap), expected)
+    assert_allclose(pymsis.calculate(date, lon, lat, alt, aps=ap), expected)
 
 
 @pytest.mark.parametrize(
@@ -401,7 +408,7 @@ def test_run_auto_f107(input_auto_f107_ap):
 def test_bad_run_inputs(inputs):
     # Inputs that have nan's in them at various places should all raise for the user
     with pytest.raises(ValueError, match="Input data has non-finite values"):
-        msis.run(*inputs)
+        pymsis.calculate(*inputs)
 
 
 @pytest.mark.parametrize(
@@ -418,7 +425,7 @@ def test_keyword_argument_call(input_data, version, msis_lib):
     dseconds = (date.astype("datetime64[s]") - date.astype("datetime64[D]")).astype(
         float
     )
-    run_output = msis.run(
+    run_output = pymsis.calculate(
         dates=date,
         alts=alt,
         lats=lat,
@@ -450,10 +457,12 @@ def test_changing_options(input_data, expected_output, expected_output_with_opti
     # so we need to make sure that we are actually changing the model
     # when the options change.
     assert_allclose(
-        np.squeeze(msis.run(*input_data, options=[1] * 25)), expected_output, rtol=1e-5
+        np.squeeze(pymsis.calculate(*input_data, options=[1] * 25)),
+        expected_output,
+        rtol=1e-5,
     )
     assert_allclose(
-        np.squeeze(msis.run(*input_data, options=[0] * 25)),
+        np.squeeze(pymsis.calculate(*input_data, options=[0] * 25)),
         expected_output_with_options,
         rtol=1e-5,
     )
@@ -468,9 +477,9 @@ def test_options_calls(input_data, version, msis_lib):
     # Reset the cache
     msis_lib._last_used_options = None
     with patch(msis_lib.__name__ + ".pyinitswitch") as mock_init:
-        msis.run(*input_data, options=[0] * 25, version=version)
+        pymsis.calculate(*input_data, options=[0] * 25, version=version)
         mock_init.assert_called_once()
-        msis.run(*input_data, options=[0] * 25, version=version)
+        pymsis.calculate(*input_data, options=[0] * 25, version=version)
         # Called again shouldn't call the initialization function
         mock_init.assert_called_once()
 
@@ -510,7 +519,7 @@ def test_multithreaded(
     def run_function(input_items):
         version, options, expected = input_items
         return np.squeeze(
-            msis.run(*input_data, version=version, options=options)
+            pymsis.calculate(*input_data, version=version, options=options)
         ), expected
 
     with concurrent.futures.ThreadPoolExecutor() as executor:
@@ -524,8 +533,8 @@ def run_function(input_items):
 
 def test_output_enum(input_data):
     # Make sure we can access the output enums
-    assert msis.Variable.MASS_DENSITY == 0
-    assert msis.Variable._member_names_ == [
+    assert pymsis.Variable.MASS_DENSITY == 0
+    assert pymsis.Variable._member_names_ == [
         "MASS_DENSITY",
         "N2",
         "O2",
@@ -538,5 +547,11 @@ def test_output_enum(input_data):
         "NO",
         "TEMPERATURE",
     ]
-    data = msis.run(*input_data)
-    assert data[..., msis.Variable.MASS_DENSITY] == data[..., 0]
+    data = pymsis.calculate(*input_data)
+    assert data[..., pymsis.Variable.MASS_DENSITY] == data[..., 0]
+
+
+def test_deprecated_legacy_imports():
+    # Make sure that msis.run() is still available and
+    # we are getting our expected variables
+    assert pymsis.calculate == msis.run == msis.calculate
diff --git a/tests/test_regression.py b/tests/test_regression.py
index 662d662..b76c494 100644
--- a/tests/test_regression.py
+++ b/tests/test_regression.py
@@ -3,7 +3,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
-from pymsis import msis
+import pymsis
 
 
 def run_input_line(line, version):
@@ -28,7 +28,7 @@ def run_input_line(line, version):
     )
 
     test_inp = (date, lon, lat, alt, f107, f107a, ap)
-    test_output = np.squeeze(msis.run(*test_inp, version=version))
+    test_output = np.squeeze(pymsis.calculate(*test_inp, version=version))
 
     # Rearrange the run's output to match the expected
     # ordering of elements