diff --git a/docs/source/reference/shared_tools/index.rst b/docs/source/reference/shared_tools/index.rst
index 053cfe34..6d232a68 100644
--- a/docs/source/reference/shared_tools/index.rst
+++ b/docs/source/reference/shared_tools/index.rst
@@ -4,12 +4,11 @@
 shared_tools
 *********************************
 
-The tools are defined in ``pyDeltaRCM.shared_tools``. 
+The tools are defined in ``pyDeltaRCM.shared_tools``.
 
+.. currentmodule:: pyDeltaRCM
 
-.. currentmodule:: pyDeltaRCM.shared_tools
+.. autosummary::
+  :toctree: ../../_autosummary
 
-.. autosummary:: 
-	:toctree: ../../_autosummary
-
-	shared_tools
+  shared_tools
diff --git a/pyDeltaRCM/init_tools.py b/pyDeltaRCM/init_tools.py
index 8cb14a7f..93867e27 100644
--- a/pyDeltaRCM/init_tools.py
+++ b/pyDeltaRCM/init_tools.py
@@ -22,6 +22,9 @@
 import logging
 import time
 import yaml
+
+from . import shared_tools
+
 # tools for initiating deltaRCM model domain
 
 
@@ -97,7 +100,7 @@ def import_files(self):
         if self.seed is not None:
             if self.verbose >= 2:
                 print("setting random seed to %s " % str(self.seed))
-            np.random.seed(self.seed)
+            shared_tools.set_random_seed(self.seed)
 
     def set_constants(self):
 
@@ -113,17 +116,13 @@ def set_constants(self):
                               [-1, 0, 1],
                               [-sqrt05, 0, sqrt05]])
 
-        self.iwalk = np.array([[-1, 0, 1],
-                               [-1, 0, 1],
-                               [-1, 0, 1]])
+        self.iwalk = shared_tools.get_iwalk()
 
         self.jvec = np.array([[-sqrt05, -1, -sqrt05],
                               [0, 0, 0],
                               [sqrt05, 1, sqrt05]])
 
-        self.jwalk = np.array([[-1, -1, -1],
-                               [0, 0, 0],
-                               [1, 1, 1]])
+        self.jwalk = shared_tools.get_jwalk()
         self.dxn_iwalk = [1, 1, 0, -1, -1, -1, 0, 1]
         self.dxn_jwalk = [0, 1, 1, 1, 0, -1, -1, -1]
         self.dxn_dist = \
diff --git a/pyDeltaRCM/sed_tools.py b/pyDeltaRCM/sed_tools.py
index c7c47084..ad384903 100644
--- a/pyDeltaRCM/sed_tools.py
+++ b/pyDeltaRCM/sed_tools.py
@@ -18,12 +18,12 @@
 
 from netCDF4 import Dataset
 
-from .shared_tools import shared_tools
+from . import shared_tools
 
 # tools for sediment routing algorithms and deposition/erosion
 
 
-class sed_tools(shared_tools):
+class sed_tools(object):
 
     def sed_route(self):
         """route all sediment"""
@@ -183,51 +183,36 @@ def sed_parcel(self, theta_sed, sed, px, py):
             # choose next with weights
 
             it += 1
-            ind = (px, py)
 
-            depth_ind = self.pad_depth[
-                ind[0] - 1 + 1:ind[0] + 2 + 1, ind[1] - 1 + 1:ind[1] + 2 + 1]
-            cell_type_ind = self.pad_cell_type[
-                ind[0] - 1 + 1:ind[0] + 2 + 1, ind[1] - 1 + 1:ind[1] + 2 + 1]
+            stage_nbrs = self.pad_stage[px - 1 + 1:px + 2 + 1, py - 1 + 1:py + 2 + 1]
+            depth_ind = self.pad_depth[px - 1 + 1:px + 2 + 1, py - 1 + 1:py + 2 + 1]
+            cell_type_ind = self.pad_cell_type[px - 1 + 1:px + 2 + 1, py - 1 + 1:py + 2 + 1]
 
-            w1 = np.maximum(0, (self.qx[ind] * self.jvec +
-                                self.qy[ind] * self.ivec))
-            w2 = depth_ind ** theta_sed
-            weight = (w1 * w2 / self.distances)
+            weights = shared_tools.get_weight_at_cell(
+                (px, py),
+                stage_nbrs.flatten(), depth_ind.flatten(), cell_type_ind.flatten(),
+                self.stage[px, py], self.qx[px, py], self.qy[px, py],
+                self.ivec.flatten(), self.jvec.flatten(), self.distances.flatten(),
+                self.dry_depth, self.gamma, theta_sed)
 
-            weight[depth_ind <= self.dry_depth] = 0.0001
-            weight[cell_type_ind == -2] = 0
+            new_cell = shared_tools.random_pick(weights)
 
-            if ind[0] == 0:
-                weight[0, :] = 0
-
-            new_cell = self.random_pick(np.cumsum(weight.flatten()))
-
-            jstep = self.iwalk.flat[new_cell]
-            istep = self.jwalk.flat[new_cell]
-            dist = np.sqrt(istep * istep + jstep * jstep)
+            dist, istep, jstep, _ = shared_tools.get_steps(new_cell, self.iwalk.flat[:], self.jwalk.flat[:])
 
             # deposition and erosion
 
             if sed == 'sand':  # sand
-                if dist > 0:
-                    # deposition in current cell
-                    self.qs[px, py] = (self.qs[px, py] +
-                                       self.Vp_res / 2 / self.dt / self.dx)
-                px = px + istep
-                py = py + jstep
-
-                if dist > 0:
-                    # deposition in downstream cell
-                    self.qs[px, py] = (self.qs[px, py] +
-                                       self.Vp_res / 2 / self.dt / self.dx)
+
+                depoPart = self.Vp_res / 2 / self.dt / self.dx
+
+                px, py, self.qs = shared_tools.partition_sand(self.qs, depoPart, py, px, dist, istep, jstep)
 
                 self.sand_dep_ero(px, py)
 
             if sed == 'mud':  # mud
 
-                px = px + istep
-                py = py + jstep
+                px = px + jstep
+                py = py + istep
 
                 self.mud_dep_ero(px, py)
 
@@ -240,8 +225,10 @@ def sand_route(self):
         theta_sed = self.theta_sand
 
         num_starts = int(self.Np_sed * self.f_bedload)
-        start_indices = [self.random_pick_inlet(
-            self.inlet) for x in range(num_starts)]
+        inlet_weights = np.ones_like(self.inlet)
+        start_indices = [
+            self.inlet[shared_tools.random_pick(inlet_weights / sum(inlet_weights))]
+            for x in range(num_starts)]
 
         for np_sed in range(num_starts):
 
@@ -281,8 +268,10 @@ def mud_route(self):
         theta_sed = self.theta_mud
 
         num_starts = int(self.Np_sed * (1 - self.f_bedload))
-        start_indices = [self.random_pick_inlet(
-            self.inlet) for x in range(num_starts)]
+        inlet_weights = np.ones_like(self.inlet)
+        start_indices = [
+            self.inlet[shared_tools.random_pick(inlet_weights / sum(inlet_weights))]
+            for x in range(num_starts)]
 
         for np_sed in range(num_starts):
 
diff --git a/pyDeltaRCM/shared_tools.py b/pyDeltaRCM/shared_tools.py
index 34ae9af6..758d95fa 100644
--- a/pyDeltaRCM/shared_tools.py
+++ b/pyDeltaRCM/shared_tools.py
@@ -3,39 +3,208 @@
 import numpy as np
 
 import time
+from numba import njit, jit, typed
 
 # tools shared between deltaRCM water and sediment routing
 
 
-class shared_tools(object):
+def get_iwalk():
+    return np.array([[-1, 0, 1],
+                     [-1, 0, 1],
+                     [-1, 0, 1]])
 
-    def random_pick(self, probs):
-        """
-        Randomly pick a number weighted by array probabilities (len 9)
-        Return the index of the selected weight in array probs
-        Takes a numpy array that is the precalculated cumulative probability
-        around the cell flattened to 1D.
-        """
 
-        idx = probs.searchsorted(np.random.uniform(0, probs[-1]))
+def get_jwalk():
+    return np.array([[-1, -1, -1],
+                     [0, 0, 0],
+                     [1, 1, 1]])
 
-        return idx
 
-    def random_pick_inlet(self, choices, probs=None):
-        """
-        Randomly pick a number from array choices weighted by array probs
-        Values in choices are column indices
+@njit
+def set_random_seed(_seed):
+    np.random.seed(_seed)
 
-        Return a tuple of the randomly picked index for row 0
-        """
 
-        if not probs:
-            probs = np.ones(len(choices))
+@njit
+def get_random_uniform(N):
+    return np.random.uniform(0, 1)
 
-        cutoffs = np.cumsum(probs)
-        idx = cutoffs.searchsorted(np.random.uniform(0, cutoffs[-1]))
 
-        return choices[idx]
+@njit
+def partition_sand(qs, depoPart, py, px, dist, istep, jstep):
+    '''Spread sand between two cells
+    '''
+    if dist > 0:
+        # deposition in current cell
+        qs[px, py] += depoPart
+
+    px = px + jstep
+    py = py + istep
+
+    if dist > 0:
+        # deposition in downstream cell
+        qs[px, py] += depoPart
+    return px, py, qs
+
+
+@njit
+def get_steps(new_cells, iwalk, jwalk):
+    '''find the values giving the next step
+    '''
+    istep = iwalk[new_cells]
+    jstep = jwalk[new_cells]
+    dist = np.sqrt(istep * istep + jstep * jstep)
+
+    astep = dist != 0
+
+    return dist, istep, jstep, astep
+
+
+@njit
+def update_dirQfield(qfield, dist, inds, astep, dirstep):
+    '''update unit vector of water flux in x or y
+    '''
+    for i, ii in enumerate(inds):
+        if astep[i]:
+            qfield[ii] += dirstep[i] / dist[i]
+    return qfield
+
+
+@njit
+def update_absQfield(qfield, dist, inds, astep, Qp_water, dx):
+    '''Update norm of water flux vector
+    '''
+    for i, ii in enumerate(inds):
+        if astep[i]:
+            qfield[ii] += Qp_water / dx / 2
+    return qfield
+
+
+@njit
+def random_pick(prob):
+    """
+    Randomly pick a number weighted by array probabilities (len 9)
+    Return the index of the selected weight in array probs
+    Takes a numpy array that is the precalculated cumulative probability
+    around the cell flattened to 1D.
+    """
+
+    arr = np.arange(len(prob))
+    return arr[np.searchsorted(np.cumsum(prob), get_random_uniform(1))]
+
+
+@njit
+def custom_unravel(i, shape):
+    """Function to unravel indexes for 2D array
+    """
+    if i > (shape[1] * shape[0]):
+        raise IndexError("Index is out of matrix bounds")
+    x = i // shape[1]
+    y = i % shape[1]
+    return x, y
+
+
+@njit
+def custom_ravel(tup, shape):
+    """Function to ravel indexes for 2D array
+    """
+    if tup[0] > shape[0] or tup[1] > shape[1]:
+        raise IndexError("Index is out of matrix bounds")
+    x = tup[0] * shape[1]
+    y = tup[1]
+    return x + y
+
+
+@njit
+def check_for_loops(indices, inds, it, L0, loopedout, domain_shape, CTR, free_surf_flag):
+
+    looped = typed.List() # numba typed list for iteration
+    for i in np.arange(len(indices)):
+        row = indices[i, :]
+        v = len(row[row > 0]) != len(set(row[row > 0]))
+        looped.append(v)
+    travel = (0, it)
+
+    for n in range(indices.shape[0]):
+        ind = inds[n]
+        if looped[n] and (ind > 0) and (max(travel) > L0):
+            loopedout[n] += 1
+            px, py = custom_unravel(ind, domain_shape)
+            Fx = px - 1
+            Fy = py - CTR
+
+            Fw = np.sqrt(Fx**2 + Fy**2)
+
+            if Fw != 0:
+                px = px + np.round(Fx / Fw * 5.)
+                py = py + np.round(Fy / Fw * 5.)
+
+            px = max(px, L0)
+            px = int(min(domain_shape[0] - 2, px))
+
+            py = max(1, py)
+            py = int(min(domain_shape[1] - 2, py))
+
+            nind = custom_ravel((px, py), domain_shape)
+
+            inds[n] = nind
+
+            free_surf_flag[n] = -1
+
+    return inds, loopedout, free_surf_flag
+
+
+@njit
+def calculate_new_ind(indices, new_cells, iwalk, jwalk, domain_shape):
+    newbies = []
+    for p, q in zip(indices, new_cells):
+        if q != 4:
+            ind_tuple = custom_unravel(p, domain_shape)
+            new_ind = (ind_tuple[0] + jwalk[q],
+                       ind_tuple[1] + iwalk[q])
+            newbies.append(custom_ravel(new_ind, domain_shape))
+        else:
+            newbies.append(0)
+
+    return np.array(newbies)
+
+
+@njit
+def get_weight_at_cell(ind, stage_nbrs, depth_nbrs, ct_nbrs, stage, qx, qy,
+                       ivec, jvec, distances, dry_depth, gamma, theta):
+
+    weight_sfc = np.maximum(0, (stage - stage_nbrs) / distances)
+
+    weight_int = np.maximum(0, (qx * jvec + qy * ivec) / distances)
+
+    if ind[0] == 0:
+        weight_sfc[:3] = np.nan
+        weight_int[:3] = np.nan
+
+    drywall = (depth_nbrs <= dry_depth) | (ct_nbrs == -2)
+    weight_sfc[drywall] = np.nan
+    weight_int[drywall] = np.nan
+
+    if np.nansum(weight_sfc) > 0:
+        weight_sfc = weight_sfc / np.nansum(weight_sfc)
+
+    if np.nansum(weight_int) > 0:
+        weight_int = weight_int / np.nansum(weight_int)
+
+    weight = gamma * weight_sfc + (1 - gamma) * weight_int
+    weight = depth_nbrs ** theta * weight
+    weight[depth_nbrs <= dry_depth] = 0
+
+    nanWeight = np.isnan(weight)
+
+    if np.any(weight[~nanWeight] != 0):
+        weight = weight / np.nansum(weight)
+        weight[nanWeight] = 0
+    else:
+        weight[~nanWeight] = 1 / len(weight[~nanWeight])
+        weight[nanWeight] = 0
+    return weight
+
 
 
 def _get_version():
diff --git a/pyDeltaRCM/water_tools.py b/pyDeltaRCM/water_tools.py
index bb26642b..0add31cf 100644
--- a/pyDeltaRCM/water_tools.py
+++ b/pyDeltaRCM/water_tools.py
@@ -1,4 +1,3 @@
-
 import sys
 import os
 import re
@@ -17,13 +16,12 @@
 from scipy import ndimage
 
 from netCDF4 import Dataset
-
-from .shared_tools import shared_tools
+from . import shared_tools
 
 # tools for water routing algorithms
 
 
-class water_tools(shared_tools):
+class water_tools(object):
 
     def init_water_iteration(self):
 
@@ -45,13 +43,16 @@ def init_water_iteration(self):
     def run_water_iteration(self):
 
         iter = 0
-        start_indices = [self.random_pick_inlet(self.inlet) for x in range(self.Np_water)]
+        inlet_weights = np.ones_like(self.inlet)
+        start_indices = [
+            self.inlet[shared_tools.random_pick(inlet_weights / sum(inlet_weights))]
+            for x in range(self.Np_water)]
 
         self.qxn.flat[start_indices] += 1
         self.qwn.flat[start_indices] += self.Qp_water / self.dx / 2
 
         self.indices[:, 0] = start_indices
-        current_inds = list(start_indices)
+        current_inds = np.array(start_indices)
 
         self.looped[:] = 0
 
@@ -70,16 +71,31 @@ def run_water_iteration(self):
                          if x != (0, 0) else 4
                          for x in inds_tuple]
 
-            new_inds = self.calculate_new_indices(inds_tuple, new_cells)
+            new_cells = np.array(new_cells)
 
-            dist = self.update_steps(current_inds, new_inds, new_cells)
+            next_index = shared_tools.calculate_new_ind(
+                current_inds,
+                new_cells,
+                self.iwalk.flatten(),
+                self.jwalk.flatten(),
+                self.eta.shape)
 
-            new_inds = np.array(new_inds, dtype=np.int)
-            new_inds[np.array(dist) == 0] = 0
+            dist, istep, jstep, astep = shared_tools.get_steps(
+                new_cells,
+                self.iwalk.flat[:],
+                self.jwalk.flat[:])
 
-            self.indices[:, iter] = new_inds
+            self.update_Q(dist, current_inds, next_index, astep, jstep, istep)
 
-            current_inds = self.check_for_loops(new_inds, iter)
+            current_inds, self.looped, self.free_surf_flag = shared_tools.check_for_loops(
+                self.indices,
+                next_index,
+                iter,
+                self.L0,
+                self.looped,
+                self.eta.shape,
+                self.CTR,
+                self.free_surf_flag)
 
             current_inds = self.check_for_boundary(current_inds)
 
@@ -193,144 +209,40 @@ def get_water_weight_array(self):
 
         for i in range(self.L):
             for j in range(self.W):
-                self.water_weights[i, j] = self.get_weight_at_cell((i, j))
-
-    def get_weight_at_cell(self, ind):
-
-        stage_ind = self.pad_stage[
-            ind[0] - 1 + 1:ind[0] + 2 + 1, ind[1] - 1 + 1:ind[1] + 2 + 1]
-
-        weight_sfc = np.maximum(0,
-                                (self.stage[ind] - stage_ind) / self.distances)
-
-        weight_int = np.maximum(0, (self.qx[ind] * self.jvec +
-                                    self.qy[ind] * self.ivec) / self.distances)
-
-        if ind[0] == 0:
-            weight_sfc[0, :] = 0
-            weight_int[0, :] = 0
-
-        depth_ind = self.pad_depth[
-            ind[0] - 1 + 1:ind[0] + 2 + 1, ind[1] - 1 + 1:ind[1] + 2 + 1]
-        ct_ind = self.pad_cell_type[
-            ind[0] - 1 + 1:ind[0] + 2 + 1, ind[1] - 1 + 1:ind[1] + 2 + 1]
-
-        weight_sfc[(depth_ind <= self.dry_depth) | (ct_ind == -2)] = 0
-        weight_int[(depth_ind <= self.dry_depth) | (ct_ind == -2)] = 0
-
-        if np.nansum(weight_sfc) > 0:
-            weight_sfc = weight_sfc / np.nansum(weight_sfc)
-
-        if np.nansum(weight_int) > 0:
-            weight_int = weight_int / np.nansum(weight_int)
-
-        weight = self.gamma * weight_sfc + (1 - self.gamma) * weight_int
-        weight = depth_ind ** self.theta_water * weight
-        weight[depth_ind <= self.dry_depth] = 0
-        if np.sum(weight) != 0:
-            weight = weight / np.sum(weight)
-        else:
-            weight = weight
-
-        return np.cumsum(weight.flatten())
+                stage_nbrs = self.pad_stage[i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]
+                depth_nbrs = self.pad_depth[i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]
+                ct_nbrs = self.pad_cell_type[i - 1 + 1:i + 2 + 1, j - 1 + 1:j + 2 + 1]
+                self.water_weights[i, j] = shared_tools.get_weight_at_cell(
+                    (i, j),
+                    stage_nbrs.flatten(), depth_nbrs.flatten(), ct_nbrs.flatten(),
+                    self.stage[i, j], self.qx[i, j], self.qy[i, j],
+                    self.ivec.flatten(), self.jvec.flatten(), self.distances.flatten(),
+                    self.dry_depth, self.gamma, self.theta_water)
 
     def get_new_cell(self, ind):
         weight = self.water_weights[ind[0], ind[1]]
-        new_cell = self.random_pick(weight)
-        return new_cell
-
-    def calculate_new_ind(self, ind, new_cell):
-
-        new_ind = (ind[0] + self.jwalk.flat[new_cell],
-                   ind[1] + self.iwalk.flat[new_cell])
-        # added wrap mode to fct to resolve ValueError due to negative numbers
-        new_ind_flat = np.ravel_multi_index(
-            new_ind, self.depth.shape, mode='wrap')
-
-        return new_ind_flat
-
-    def calculate_new_indices(self, inds_tuple, new_cells):
-        newbies = []
-        for i, j in zip(inds_tuple, new_cells):
-            if j != 4:
-                newbies.append(self.calculate_new_ind(i, j))
-            else:
-                newbies.append(0)
-        return newbies
-
-    def step_update(self, ind, new_ind, new_cell):
-
-        istep = self.iwalk.flat[new_cell]
-        jstep = self.jwalk.flat[new_cell]
-        dist = np.sqrt(istep * istep + jstep * jstep)
-
-        if dist > 0:
-
-            self.qxn.flat[ind] += jstep / dist
-            self.qyn.flat[ind] += istep / dist
-            self.qwn.flat[ind] += self.Qp_water / self.dx / 2
-
-            self.qxn.flat[new_ind] += jstep / dist
-            self.qyn.flat[new_ind] += istep / dist
-            self.qwn.flat[new_ind] += self.Qp_water / self.dx / 2
-
-        return dist
-
-    def update_steps(self, current_inds, new_inds, new_cells):
-        newbies = []
-        for i, j, k in zip(current_inds, new_inds, new_cells):
-            if i > 0:
-                newbies.append(self.step_update(i, j, k))
-            else:
-                newbies.append(0)
-        return newbies
-
-    def check_for_loops(self, inds, it):
-
-        looped = [len(i[i > 0]) != len(set(i[i > 0])) for i in self.indices]
-
-        for n in range(self.Np_water):
-
-            ind = inds[n]
-            # revised 'it' to 'np.max(it)' to match python 2 > assessment
-            if (looped[n]) and (ind > 0) and (np.max(it) > self.L0):
-
-                self.looped[n] += 1
-
-                it = np.unravel_index(ind, self.depth.shape)
-
-                px, py = it
-
-                Fx = px - 1
-                Fy = py - self.CTR
-
-                Fw = np.sqrt(Fx**2 + Fy**2)
-
-                if Fw != 0:
-                    px = px + np.round(Fx / Fw * 5.)
-                    py = py + np.round(Fy / Fw * 5.)
-
-                px = max(px, self.L0)
-                px = int(min(self.L - 2, px))
-
-                py = max(1, py)
-                py = int(min(self.W - 2, py))
-
-                nind = np.ravel_multi_index((px, py), self.depth.shape)
-
-                inds[n] = nind
-
-                self.free_surf_flag[n] = -1
-
-        return inds
+        return shared_tools.random_pick(weight)
+
+    def update_Q(self, dist, current_inds, next_index, astep, jstep, istep):
+
+        self.qxn = shared_tools.update_dirQfield(self.qxn.flat[:], dist, current_inds, astep, jstep
+                                          ).reshape(self.qxn.shape)
+        self.qyn = shared_tools.update_dirQfield(self.qyn.flat[:], dist, current_inds, astep, istep
+                                          ).reshape(self.qyn.shape)
+        self.qwn = shared_tools.update_absQfield(self.qwn.flat[:], dist, current_inds, astep, self.Qp_water, self.dx
+                                          ).reshape(self.qwn.shape)
+        self.qxn = shared_tools.update_dirQfield(self.qxn.flat[:], dist, next_index, astep, jstep
+                                          ).reshape(self.qxn.shape)
+        self.qyn = shared_tools.update_dirQfield(self.qyn.flat[:], dist, next_index, astep, istep
+                                          ).reshape(self.qyn.shape)
+        self.qwn = shared_tools.update_absQfield(self.qwn.flat[:], dist, next_index, astep, self.Qp_water, self.dx
+                                          ).reshape(self.qwn.shape)
 
     def check_for_boundary(self, inds):
 
-        self.free_surf_flag[(self.cell_type.flat[inds] == -1)
-                            & (self.free_surf_flag == 0)] = 1
+        self.free_surf_flag[(self.cell_type.flat[inds] == -1) & (self.free_surf_flag == 0)] = 1
 
-        self.free_surf_flag[(self.cell_type.flat[inds] == -1)
-                            & (self.free_surf_flag == -1)] = 2
+        self.free_surf_flag[(self.cell_type.flat[inds] == -1) & (self.free_surf_flag == -1)] = 2
 
         inds[self.free_surf_flag == 2] = 0
 
diff --git a/requirements.txt b/requirements.txt
index 90684434..d60e3ff1 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,3 +4,4 @@ scipy
 netCDF4
 basic_modeling_interface
 pyyaml
+numba
diff --git a/setup.py b/setup.py
index 3bb8c793..97f73513 100644
--- a/setup.py
+++ b/setup.py
@@ -13,7 +13,8 @@
       include_package_data=True,
       url='https://github.com/DeltaRCM/pyDeltaRCM',
       install_requires=['matplotlib', 'netCDF4',
-                        'basic-modeling-interface', 'scipy', 'numpy', 'pyyaml'],
+                        'basic-modeling-interface', 'scipy', 'numpy', 'pyyaml',
+                        'numba'],
       entry_points={
                   'console_scripts': ['run_pyDeltaRCM=pyDeltaRCM.command_line:run_model'],
       }
diff --git a/tests/test_deltaRCM_driver.py b/tests/test_deltaRCM_driver.py
index ea7c2408..92440731 100644
--- a/tests/test_deltaRCM_driver.py
+++ b/tests/test_deltaRCM_driver.py
@@ -20,17 +20,15 @@ def test_init():
     assert delta._is_finalized == False
 
 
-delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
-
-
 def test_update():
-    
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
     delta.update()
     assert delta._time == 1.0
     assert delta._is_finalized == False
 
 
 def test_finalize():
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
     delta.finalize()
     assert delta._is_finalized == True
 
@@ -40,9 +38,60 @@ def test_multifinalization_error():
     err_delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
     err_delta.update()
     # test will Fail if any assertion is wrong
-    assert err_delta._time == 1.0 
+    assert err_delta._time == 1.0
     assert err_delta._is_finalized == False
     err_delta.finalize()
     assert err_delta._is_finalized == True
     # next line should throw RuntimeError and test will xFail
     err_delta.update()
+
+
+def test_getters_setters():
+
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+
+    assert np.all(delta.sea_surface_elevation == 0)
+    assert delta.water_depth[0, 2] == 0
+    assert delta.water_depth[0, 3] == 1
+    assert delta.water_depth[4, 4] == 1
+    assert delta.bed_elevation[0, 2] == 0
+    assert delta.bed_elevation[0, 3] == -1
+    assert delta.bed_elevation[4, 4] == -1
+
+    assert delta.sea_surface_mean_elevation == 0
+    delta.sea_surface_mean_elevation = 0.5
+    assert delta.sea_surface_mean_elevation == 0.5
+
+    assert delta.sea_surface_elevation_change == 0.001
+    delta.sea_surface_elevation_change = 0.002
+    assert delta.sea_surface_elevation_change == 0.002
+
+    assert delta.bedload_fraction == 0.5
+    delta.bedload_fraction = 0.25
+    assert delta.bedload_fraction == 0.25
+
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+
+    assert delta.channel_flow_velocity == 1
+    delta.channel_flow_velocity = 3
+    assert delta.channel_flow_velocity == 3
+    assert delta.u_max == 6
+
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+
+    assert delta.channel_width == 2
+    delta.channel_width = 10
+    assert delta.channel_width == 10
+    assert delta.N0 == 10
+
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+
+    assert delta.channel_flow_depth == 1
+    delta.channel_flow_depth = 2
+    assert delta.channel_flow_depth == 2
+
+    delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
+
+    assert delta.influx_sediment_concentration == 0.1
+    delta.influx_sediment_concentration = 2
+    assert delta.C0 == 0.02
diff --git a/tests/test_sed_tools.py b/tests/test_sed_tools.py
index fb577328..29c21e06 100644
--- a/tests/test_sed_tools.py
+++ b/tests/test_sed_tools.py
@@ -21,6 +21,8 @@ def test_sed_route():
     test the function sed_tools.sed_route
     """
     delta.pad_cell_type = np.pad(delta.cell_type, 1, 'edge')
+    delta.pad_stage = np.pad(delta.stage, 1, 'edge')
+    delta.pad_depth = np.pad(delta.depth, 1, 'edge')
     delta.sed_route()
     [a, b] = np.shape(delta.pad_depth)
 
diff --git a/tests/test_shared_tools.py b/tests/test_shared_tools.py
index 6580cd57..d29d496c 100644
--- a/tests/test_shared_tools.py
+++ b/tests/test_shared_tools.py
@@ -16,6 +16,85 @@
 # delta._delta.**shared_tools_function**
 
 
+def test_set_random_assignments():
+    """
+    Test for function shared_tools.get_random_uniform and
+    test for function shared_tools.set_random_seed
+    """
+    shared_tools.set_random_seed(delta.seed)
+    got = shared_tools.get_random_uniform(1)
+    _exp = 0.5488135039273248
+    assert got == pytest.approx(_exp)
+
+
+def test_sand_partition():
+    """
+    Test for function shared_tools.partition_sand
+    """
+    nx, ny, qsn = shared_tools.partition_sand(
+        delta.qs, 1, 4, 4, 1, 0, 1
+    )
+    assert nx == 5
+    assert ny == 4
+    assert qsn[4, 4] == 1
+    assert qsn[5, 4] == 1
+    assert np.all(qsn[qsn != 1] == 0)
+
+
+def test_get_steps():
+    """
+    Test for function shared_tools.get_steps
+    """
+    new_cells = np.arange(9)
+    iwalk = shared_tools.get_iwalk()
+    jwalk = shared_tools.get_jwalk()
+
+    d, i, j, a = shared_tools.get_steps(new_cells, iwalk.flatten(), jwalk.flatten())
+
+    d_exp = np.array([1.41421356, 1., 1.41421356, 1., 0.,
+             1., 1.41421356, 1., 1.41421356])
+    i_exp = np.array([-1,  0,  1, -1,  0,  1, -1,  0,  1])
+    j_exp = np.array([-1, -1, -1,  0,  0,  0,  1,  1,  1])
+
+    assert np.all(np.delete(a, 4))
+    assert ~a[4]
+    assert i == pytest.approx(i_exp)
+    assert j == pytest.approx(j_exp)
+    assert d == pytest.approx(d_exp)
+
+
+def test_update_dirQfield():
+    """
+    Test for function shared_tools.update_dirQfield
+    """
+    np.random.seed(delta.seed)
+    qx = np.random.uniform(0, 10, 9)
+    d = np.array([1, np.sqrt(2), 0])
+    astep = np.array([True, True, False])
+    inds = np.array([3, 4, 5])
+    stepdir = np.array([1, 1, 0])
+    qxn = shared_tools.update_dirQfield(np.copy(qx), d, inds, astep, stepdir)
+    qxdiff = qxn - qx
+    qxdiff_exp = np.array([1, np.sqrt(2) / 2, 0])
+    assert np.all(qxdiff[3:6] == pytest.approx(qxdiff_exp))
+
+
+def test_update_absQfield():
+    """
+    Test for function shared_tools.update_absQfield
+    """
+    np.random.seed(delta.seed)
+    qw = np.random.uniform(0, 10, 9)
+    d = np.array([1, np.sqrt(2), 0])
+    astep = np.array([True, True, False])
+    inds = np.array([3, 4, 5])
+    qwn = shared_tools.update_absQfield(np.copy(qw), d, inds, astep, delta.Qp_water, delta.dx)
+    qwdiff = qwn - qw
+    diffelem = delta.Qp_water / delta.dx / 2
+    qwdiff_exp = np.array([diffelem, diffelem, 0])
+    assert np.all(qwdiff[3:6] == pytest.approx(qwdiff_exp))
+
+
 def test_random_pick():
     """
     Test for function shared_tools.random_pick
@@ -24,13 +103,129 @@ def test_random_pick():
     probs = np.zeros((8,))
     probs[0] = 1
     # should return first index
-    assert delta.random_pick(probs) == 0
+    assert shared_tools.random_pick(probs) == 0
 
 
-def test_random_pick_inlet():
+def test_random_pick_anybut_first():
     """
-    Test for function shared_tools.random_pick_inlet
+    Test for function shared_tools.random_pick
     """
-    choices = [0]
-    probs = np.ones((1,))
-    assert delta.random_pick_inlet(choices, probs) == 0
+    shared_tools.set_random_seed(delta.seed)
+    probs = (1 / 7) * np.ones((3, 3), dtype=np.float64)
+    probs[0, 0] = 0
+    probs[1, 1] = 0
+    # should never return first index
+    _rets = np.zeros((100,))
+    for i in range(100):
+        _rets[i] = shared_tools.random_pick(probs.flatten())
+    assert np.all(_rets != 0)
+    assert np.all(_rets != 4)  # THIS LINE NEEDS TO PASS!!
+    assert np.sum(_rets == 1) > 0
+    assert np.sum(_rets == 2) > 0
+    assert np.sum(_rets == 3) > 0
+    assert np.sum(_rets == 5) > 0
+    assert np.sum(_rets == 6) > 0
+    assert np.sum(_rets == 7) > 0
+    assert np.sum(_rets == 8) > 0
+
+
+def test_custom_unravel_square():
+    arr = np.arange(9).reshape((3, 3))
+    # test upper left corner
+    x, y = shared_tools.custom_unravel(0, arr.shape)
+    assert x == 0
+    assert y == 0
+    # test center
+    x, y = shared_tools.custom_unravel(4, arr.shape)
+    assert x == 1
+    assert y == 1
+    # test off-center
+    x, y = shared_tools.custom_unravel(5, arr.shape)
+    assert x == 1
+    assert y == 2
+    # test lower right corner
+    x, y = shared_tools.custom_unravel(8, arr.shape)
+    assert x == 2
+    assert y == 2
+
+
+def test_custom_unravel_rectangle():
+    arr = np.arange(50).reshape((5, 10))
+    # test a few spots
+    x, y = shared_tools.custom_unravel(19, arr.shape)
+    assert x == 1
+    assert y == 9
+    x, y = shared_tools.custom_unravel(34, arr.shape)
+    assert x == 3
+    assert y == 4
+
+
+@pytest.mark.xfail(raises=IndexError, strict=True)
+def test_custom_unravel_exceed_error():
+    arr = np.arange(9).reshape((3, 3))
+    x, y = shared_tools.custom_unravel(99, arr.shape)
+
+
+def test_check_for_loops():
+
+    idxs = np.array(
+        [[0, 11, 12, 13, 23, 22, 12],
+         [0, 1, 2, 3, 4, 5, 16]])
+    nidx = np.array([21, 6])
+    itt = 6
+    free = np.array([1, 1])
+    CTR = 4
+    L0 = 1
+    looped = np.array([0, 0])
+
+    nidx, looped, free = shared_tools.check_for_loops(
+        idxs, nidx, itt, L0, looped, (10, 10), CTR, free)
+
+    assert np.all(nidx == [41, 6])
+    assert np.all(looped == [1, 0])
+    assert np.all(free == [-1, 1])
+
+
+def test_calculate_new_ind():
+
+    cidx = np.array([12, 16, 16])
+    ncel = np.array([6, 1, 4])
+    iwalk = shared_tools.get_iwalk()
+    jwalk = shared_tools.get_jwalk()
+
+    nidx = shared_tools.calculate_new_ind(cidx, ncel, iwalk.flatten(), jwalk.flatten(), (10, 10))
+
+    nidx_exp = np.array([21, 6, 0])
+    assert np.all(nidx == nidx_exp)
+
+
+def test_get_weight_at_cell():
+
+    ind = (0, 4)
+    np.random.seed(delta.seed)
+    stage = np.random.uniform(0.5, 1, 9)
+    eta = np.random.uniform(0, 0.85, 9)
+    depth = stage - eta
+    depth[depth < 0] = 0
+    celltype = np.array([-2, -2, -2, 1, 1, -2, 0, 0, 0])
+    qx = 1
+    qy = 1
+    ivec = delta.ivec.flatten()
+    jvec = delta.jvec.flatten()
+    dists = delta.distances.flatten()
+    dry_thresh = 0.1
+    gamma = 0.001962
+    theta = 1
+
+    wts = shared_tools.get_weight_at_cell(ind, stage, depth, celltype, stage, qx, qy,
+                           ivec, jvec, dists, dry_thresh, gamma, theta)
+    assert np.all(wts[[0, 1, 2, 5]] == 0)
+    assert wts[4] == 0
+    assert np.any(wts[[3, 6, 7, 8]] != 0)
+
+
+def test_version_is_valid():
+    v = shared_tools._get_version()
+    assert type(v) is str
+    dots = [i for i, c in enumerate(v) if c == '.']
+    assert len(dots) == 2
diff --git a/tests/test_water_tools.py b/tests/test_water_tools.py
index c6643847..e95bbd0d 100644
--- a/tests/test_water_tools.py
+++ b/tests/test_water_tools.py
@@ -8,6 +8,7 @@
 
 from pyDeltaRCM.deltaRCM_driver import pyDeltaRCM
 from pyDeltaRCM import water_tools
+from pyDeltaRCM import shared_tools
 
 # need to create a simple case of pydeltarcm object to test these functions
 delta = pyDeltaRCM(input_file=os.path.join(os.getcwd(), 'tests', 'test.yaml'))
@@ -76,8 +77,9 @@ def test_calculate_new_ind():
     Test for function water_tools.calculate_new_ind
     """
     # assign old index
-    old_ind = [0, 4]
+    old_inds = np.array([4, 5])
     # assign new cell location
-    new_cell = 7
+    new_cells = np.array([7, 7])
     # expect new cell to be in location (1,4) -> 14
-    assert delta.calculate_new_ind(old_ind, new_cell) == 14
+    new_inds = shared_tools.calculate_new_ind(old_inds, new_cells, delta.iwalk.flatten(), delta.jwalk.flatten(), delta.eta.shape)
+    assert np.all(new_inds == np.array([14, 15]))
diff --git a/tests/test_yaml_parsing.py b/tests/test_yaml_parsing.py
index 0703ec29..90da48ce 100644
--- a/tests/test_yaml_parsing.py
+++ b/tests/test_yaml_parsing.py
@@ -6,6 +6,7 @@
 
 from pyDeltaRCM.deltaRCM_driver import pyDeltaRCM
 
+from pyDeltaRCM.shared_tools import set_random_seed, get_random_uniform
 
 # utilities for file writing
 def create_temporary_file(tmp_path, file_name):
@@ -105,13 +106,13 @@ def test_random_seed_settings_value(tmp_path):
     p, f = create_temporary_file(tmp_path, file_name)
     write_parameter_to_file(f, 'seed', 9999)
     f.close()
-    np.random.seed(9999)
-    _preval_same = np.random.uniform()
-    np.random.seed(5)
-    _preval_diff = np.random.uniform(1000)
+    set_random_seed(9999)
+    _preval_same = get_random_uniform(1)
+    set_random_seed(5)
+    _preval_diff = get_random_uniform(1000)
     delta = pyDeltaRCM(input_file=p)
     assert delta.seed == 9999
-    _postval_same = np.random.uniform()
+    _postval_same = get_random_uniform(1)
     assert _preval_same == _postval_same
     assert delta.seed == 9999