diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml
new file mode 100644
index 0000000..1294d6d
--- /dev/null
+++ b/.github/workflows/pypi.yml
@@ -0,0 +1,55 @@
+name: PyPI
+
+on:
+  pull_request:
+  push:
+    branches:
+      - main
+  release:
+    types:
+      - published
+
+defaults:
+  run:
+    shell: bash
+
+jobs:
+  packages:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v3
+      with:
+      # Should be enough for setuptools-scm
+      fetch-depth: 100
+      persist-credentials: false
+
+    - name: Get tags
+      run: git fetch origin 'refs/tags/*:refs/tags/*'
+
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: 3.x
+
+    - name: Install build tools
+      run: |
+        python -m pip install --upgrade pip build
+
+    - name: Build binary wheel
+      run: python -m build --sdist --wheel . --outdir dist
+
+    - name: CheckFiles
+      run: |
+        ls -lh dist
+
+    - name: Test wheels
+      run: |
+        cd dist && python -m pip install odvc*.whl
+        python -m twine check *
+
+    - name: Publish a Python distribution to PyPI
+      if: success() && github.event_name == 'release'
+      uses: pypa/gh-action-pypi-publish@master
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_PASSWORD }}
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b89eecc
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+_version.py
+*.egg-info/
+.coverage
+.ipynb_checkpoints/
+build/
+dist/
diff --git a/odvc/__init__.py b/odvc/__init__.py
index 90f4ddb..9dfa33d 100644
--- a/odvc/__init__.py
+++ b/odvc/__init__.py
@@ -1,3 +1,5 @@
+"""Ocean Dimensionless Vertical Coordinates."""
+
 try:
     from ._version import __version__
 except ImportError:
diff --git a/odvc/formulas.py b/odvc/formulas.py
index 1854cdf..d55d0ce 100644
--- a/odvc/formulas.py
+++ b/odvc/formulas.py
@@ -1,9 +1,11 @@
+"""Ocean Dimensionless Vertical Coordinates formulas."""
+
 import numpy as np
 
 
 def ocean_double_sigma_coordinate(sigma, depth, z1, z2, a, href, k_c):
     """
-    Creates a dimensioned version of ocean double sigma.
+    Create a dimensioned version of ocean double sigma.
 
     Definition:
       for k <= k_c
@@ -24,7 +26,7 @@ def ocean_double_sigma_coordinate(sigma, depth, z1, z2, a, href, k_c):
 
 def ocean_sigma_z_coordinate(sigma, eta, depth, depth_c, nsigma, zlev):
     """
-    Creates a dimensioned version of ocean sigma over z.
+    Create a dimensioned version of ocean sigma over z.
 
     Definition:
         for k <= nsigma
@@ -46,7 +48,7 @@ def ocean_sigma_z_coordinate(sigma, eta, depth, depth_c, nsigma, zlev):
 
 def ocean_sigma_coordinate(sigma, eta, depth):
     """
-    Creates a dimensioned version of ocean sigma coordinate.
+    Create a dimensioned version of ocean sigma coordinate.
 
     z(n, k, j, i) = eta(n, j, i) + sigma(k) *
                     (depth(j, i) + eta(n, j, i))
@@ -57,7 +59,7 @@ def ocean_sigma_coordinate(sigma, eta, depth):
 
 def ocean_s_coordinate(s, eta, depth, a, b, depth_c):
     """
-    Creates a dimensioned version of ocean s-coordinate.
+    Create a dimensioned version of ocean s-coordinate.
 
     z(n,k,j,i) = eta(n,j,i)*(1+s(k)) + depth_c*s(k) +
                     (depth(j,i)-depth_c)*C(k)
@@ -75,7 +77,7 @@ def ocean_s_coordinate(s, eta, depth, a, b, depth_c):
 
 def ocean_s_coordinate_g1(s, c, eta, depth, depth_c):
     """
-    Creates a dimensioned version of ocean s-coordinate generic form 1.
+    Create a dimensioned version of ocean s-coordinate generic form 1.
 
     z(n,k,j,i) = S(k,j,i) + eta(n,j,i) * (1 + S(k,j,i) / depth(j,i))
 
@@ -89,7 +91,7 @@ def ocean_s_coordinate_g1(s, c, eta, depth, depth_c):
 
 def ocean_s_coordinate_g2(s, eta, depth, depth_c, c):
     """
-    Creates a dimensioned version of s-coordinate generic form 2.
+    Create a dimensioned version of s-coordinate generic form 2.
 
     z(n,k,j,i) = eta(n,j,i) + (eta(n,j,i) + depth(j,i)) * S(k,j,i)
 
diff --git a/odvc/utils.py b/odvc/parse_formula_terms.py
similarity index 80%
rename from odvc/utils.py
rename to odvc/parse_formula_terms.py
index 66a73cd..f76471a 100644
--- a/odvc/utils.py
+++ b/odvc/parse_formula_terms.py
@@ -1,17 +1,14 @@
+"""Functions to parse the formula terms and prepare arrays."""
+
 from collections import OrderedDict
 
-import dask
-import dask.array as da
 import netCDF4
 
-# Remote THREDDS access must be single threaded.
-# This should be an option in `odvc` so local files are not penalized too.
-dask.config.set(scheduler="single-threaded")
-
 
 def get_formula_terms(var):
     """
     Return a `formula_terms` dict mapping var_names to variables.
+
     The input can be a netCDF variable (`var`) holding the `formula_terms`
     attribute or the attribute itself.
 
@@ -26,11 +23,7 @@ def get_formula_terms(var):
 
 
 def get_formula_terms_variables(nc):
-    """
-    Return a list with all variables from the `nc` object that holds the
-    `formula_terms` attribute.
-
-    """
+    """Return a list with all variables that holds the formula_terms attribute."""
 
     def func(v):
         return v is not None
@@ -46,11 +39,7 @@ def func(v):
 
 
 def get_formula_terms_dims(nc, formula_terms):
-    """
-    Returns an OrderedDict object `dims` holding the `formula_terms` dimensions
-    listed in the netCDF4-python object `nc`
-
-    """
+    """Return an OrderedDict object `dims` holding the `formula_terms` dimensions."""
     dims = OrderedDict()
     for k, v in formula_terms.items():
         dims.update({k: nc[v].dimensions})
@@ -58,11 +47,7 @@ def get_formula_terms_dims(nc, formula_terms):
 
 
 def z_shape(nc, dims):
-    """
-    Returns the vertical coordinate `shape` based on
-    the combined dimensions of the formula_terms `dims`.
-
-    """
+    """Return the vertical coordinate `shape` based on the combined dimensions of the formula_terms `dims`."""
     all_dims = (
         dims.get("eta"),
         dims.get("depth"),
@@ -87,7 +72,8 @@ def z_shape(nc, dims):
     return tuple(shape)
 
 
-def reshape(arr, new_shape):
+def _reshape(arr, new_shape):
+    """Reshape arrays according to the expected dimensions."""
     shape = arr.shape
     dims = [k for k in shape]
     slicer = slice(None, None, None)
@@ -101,6 +87,9 @@ def select(dim):
 
 
 def prepare_arrays(nc, formula_terms, new_shape):
+    """Prepare arrays for formulas."""
+    import dask.array as da
+
     arrays = {}
     for term, var in formula_terms.items():
         var = nc[var]
@@ -116,20 +105,23 @@ def prepare_arrays(nc, formula_terms, new_shape):
             if term == "eta" and var.ndim == 2:
                 chunks = (1,) + var.shape[1:]
             arr = da.from_array(var, chunks=chunks)
-            arr = reshape(arr, new_shape)
+            arr = _reshape(arr, new_shape)
         arrays.update({term: arr})
     return arrays
 
 
 def _filter_none(seq):
+    """Filter out None in a list."""
     return [x for x in seq if x is not None]
 
 
 def _flatten(seq):
+    """Flatten a list."""
     return [item for sublista in seq for item in sublista]
 
 
 def _remove_duplicate(seq):
+    """Drop duplicated from a list."""
     seen = set()
     seen_add = seen.add
     return [x for x in seq if not (x in seen or seen_add(x))]
diff --git a/pyproject.toml b/pyproject.toml
index a50133b..92fc528 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,5 @@
 [build-system]
-requires = ["setuptools>=41.2", "setuptools_scm", "wheel"]
+requires = ["setuptools>=41.2", "setuptools_scm", "wheel", "twine"]
 build-backend = "setuptools.build_meta"
 
 [tool.interrogate]
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 8512c65..2eea4f0 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -1,9 +1,11 @@
 black
+dask
 flake8-builtins
 flake8-comprehensions
 flake8-mutable
 flake8-print
 isort
+matplotlib
 nbsphinx
 pycodestyle
 pylint
diff --git a/setup.cfg b/setup.cfg
index 5498fc4..7a5c164 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,7 +30,6 @@ classifiers =
 zip_safe = True
 include_package_data = True
 install_requires =
-    dask >=1
     netcdf4
     numpy
 python_requires = >=3.7
@@ -43,7 +42,6 @@ formats = gztar
 ignore =
     *.yml
     .coveragerc
-    Makefile
     docs
     docs/*
     notebooks
diff --git a/tests/test_ocean_s_coordinate_g1.py b/tests/test_ocean_s_coordinate_g1.py
index 76819c4..3137d3e 100644
--- a/tests/test_ocean_s_coordinate_g1.py
+++ b/tests/test_ocean_s_coordinate_g1.py
@@ -1,3 +1,5 @@
+"""Test ocean_s_coordinate_g1"""
+
 import os
 import unittest
 
@@ -6,7 +8,7 @@
 from netCDF4 import Dataset
 
 from odvc import ocean_s_coordinate_g1
-from odvc.utils import (
+from odvc.parse_formula_terms import (
     get_formula_terms,
     get_formula_terms_dims,
     get_formula_terms_variables,
diff --git a/tests/test_utils.py b/tests/test_parse_formula_terms.py
similarity index 98%
rename from tests/test_utils.py
rename to tests/test_parse_formula_terms.py
index c520b27..4fa6b90 100644
--- a/tests/test_utils.py
+++ b/tests/test_parse_formula_terms.py
@@ -6,7 +6,7 @@
 import numpy as np
 import pytest
 
-from odvc.utils import (
+from odvc.parse_formula_terms import (
     get_formula_terms,
     get_formula_terms_dims,
     get_formula_terms_variables,