diff --git a/conda_package/docs/api.rst b/conda_package/docs/api.rst
index 95433fd0f..a14f213d5 100644
--- a/conda_package/docs/api.rst
+++ b/conda_package/docs/api.rst
@@ -170,6 +170,19 @@ CIME constants
 
    constants
 
+Landice Tools
+=============
+
+.. currentmodule:: mpas_tools.landice
+
+.. autosummary::
+   :toctree: generated/
+
+    visualization
+    visualization.plot_transect
+    visualization.plot_map
+    visualization.plot_grounding_lines
+
 Ocean Tools
 ===========
 
diff --git a/conda_package/docs/images/ais_map.png b/conda_package/docs/images/ais_map.png
new file mode 100644
index 000000000..e3539fb07
Binary files /dev/null and b/conda_package/docs/images/ais_map.png differ
diff --git a/conda_package/docs/images/ais_map_with_grounding_lines.png b/conda_package/docs/images/ais_map_with_grounding_lines.png
new file mode 100644
index 000000000..1df8b5c51
Binary files /dev/null and b/conda_package/docs/images/ais_map_with_grounding_lines.png differ
diff --git a/conda_package/docs/images/ais_thickness_change_map.png b/conda_package/docs/images/ais_thickness_change_map.png
new file mode 100644
index 000000000..18378d90a
Binary files /dev/null and b/conda_package/docs/images/ais_thickness_change_map.png differ
diff --git a/conda_package/docs/images/thwaites_temperature_transects.png b/conda_package/docs/images/thwaites_temperature_transects.png
new file mode 100644
index 000000000..1c28f4b0c
Binary files /dev/null and b/conda_package/docs/images/thwaites_temperature_transects.png differ
diff --git a/conda_package/docs/images/thwaites_transect.png b/conda_package/docs/images/thwaites_transect.png
new file mode 100644
index 000000000..a44bae0d4
Binary files /dev/null and b/conda_package/docs/images/thwaites_transect.png differ
diff --git a/conda_package/docs/index.rst b/conda_package/docs/index.rst
index f2866cd66..b2deb349b 100644
--- a/conda_package/docs/index.rst
+++ b/conda_package/docs/index.rst
@@ -35,6 +35,12 @@ analyzing simulations, and in other MPAS-related workflows.
 
    visualization
 
+.. toctree::
+   :caption: Landice Tools
+   :maxdepth: 2
+
+   landice/visualization
+
 .. toctree::
    :caption: Ocean Tools
    :maxdepth: 2
diff --git a/conda_package/docs/landice/visualization.rst b/conda_package/docs/landice/visualization.rst
new file mode 100644
index 000000000..7cffd37dd
--- /dev/null
+++ b/conda_package/docs/landice/visualization.rst
@@ -0,0 +1,121 @@
+.. _landice_visualization:
+
+*************
+Visualization
+*************
+
+.. _landice_transects:
+
+Plotting Landice Transects
+==========================
+
+.. image:: ../images/thwaites_transect.png
+   :width: 500 px
+   :align: center
+
+The function :py:func:`mpas_tools.landice.visualization.plot_transect()` can be
+used to plot transects of various cell-centered MALI variables. The above figure
+was plotted with the code:
+
+.. code-block:: python
+
+    import matplotlib.pyplot as plt
+    from netCDF4 import Dataset
+    from mpas_tools.landice.visualization import plot_transect
+
+
+    fig, ax = plt.subplots(2,1, figsize=(5,5), sharex=True)
+    data_path = "ISMIP6/archive/expAE02/"
+    file = "2200.nc"
+    x = [-1589311.171,-1208886.353]
+    y = [-473916.7182,-357411.6176]
+    times = [0, 1, 2, 3, 4]
+    plot_transect(data_path + file, "geometry", ax=ax[0], times=times,  x=x, y=y)
+    plot_transect(data_path + file, "surfaceSpeed", ax=ax[1], times=times,  x=x, y=y)
+    fig.savefig("thwaites_transect.png", dpi=400)
+
+This can also be used to produce transects of ice temperature, with one
+time level should be plotted per subplot:
+
+.. code-block:: python
+
+    fig, axs = plt.subplots(2,2, figsize=(10,10), sharex=True, sharey=True)
+    files = ["rst.2015-01-01.nc", "rst.2100-01-01.nc", "rst.2200-01-01.nc", "rst.2300-01-01.nc"]
+    x = [-1589311.171,-1208886.353]
+    y = [-473916.7182,-357411.6176]
+
+    for ax, file in zip(axs.ravel(), files):
+        plot_transect(data_path + file, "temperature", ax=ax,
+                      times=[0], x=x, y=y)
+        ax.set_title(file[4:8])
+    fig.savefig("thwaites_temperature_transects.png", dpi=400)
+
+.. image:: ../images/thwaites_temperature_transects.png
+   :width: 500 px
+   :align: center
+
+Plotting Landice Maps
+=====================
+
+.. image:: ../images/ais_map.png
+   :width: 500 px
+   :align: center
+
+The function :py:func:`mpas_tools.landice.visualization.plot_map()` can be
+used to plot maps of cell-centered MALI variables. The above
+figure was plotted using:
+
+.. code-block:: python
+
+    from mpas_tools.landice.visualization import plot_map
+
+
+    fig, ax = plt.subplots(1, 1, figsize=(10,8), layout='constrained')
+    var_plot, cbar, gl_plot = plot_map(data_path + files[0], "bedTopography",
+                                       ax=ax, time=0, plot_grounding_line=True)
+
+This can also be used to plot derived fields, such as thickness change from 2015 to 2300:
+
+.. code-block:: python
+
+    data_2300 = Dataset(data_path + "output_state_2300.nc")
+    data_2300.set_auto_mask(False)
+    data_2015 = Dataset(data_path + "output_state_2015.nc")
+    data_2015.set_auto_mask(False)
+    thk_2300 = data_2300.variables["thickness"][:]
+    thk_2015 = data_2015.variables["thickness"][:]
+
+    fig, ax = plt.subplots(1, 1, figsize=(6,5), layout='constrained')
+    thk_diff_plot, thk_diff_cbar, _ = plot_map(
+        data_path + "output_state_2015.nc", variable=thk_2300-thk_2015,
+        ax=ax, time=0, vmin=-1000, vmax=1000, plot_grounding_line=False,
+        variable_name="Thickness change 2015 to 2300 (m)")
+    ax.grid(linestyle='dashed')
+
+.. image:: ../images/ais_thickness_change_map.png
+   :width: 500 px
+   :align: center
+
+Plotting Landice Grounding Lines
+================================
+
+.. image:: ../images/ais_map_with_grounding_lines.png
+   :width: 500 px
+   :align: center
+
+The function :py:func:`mpas_tools.landice.visualization.plot_grounding_lines()` can be
+used to add grounding lines to maps as countours.
+
+.. code-block:: python
+
+    from mpas_tools.landice.visualization import plot_grounding_lines
+
+
+    fig, ax = plt.subplots(1, 1, figsize=(8,6), layout='constrained')
+    plot_grounding_lines([data_path + "output_state_2015.nc",
+                          data_path + "output_state_2100.nc",
+                          data_path + "output_state_2200.nc",
+                          data_path + "output_state_2300.nc"],
+                         ax=ax,
+                         cmap="plasma", times=[0])
+    plot_map(data_path + files[0], "bedTopography", ax=ax, time=0)
diff --git a/conda_package/mpas_tools/landice/__init__.py b/conda_package/mpas_tools/landice/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/conda_package/mpas_tools/landice/visualization.py b/conda_package/mpas_tools/landice/visualization.py
new file mode 100644
index 000000000..c90506841
--- /dev/null
+++ b/conda_package/mpas_tools/landice/visualization.py
@@ -0,0 +1,520 @@
+import numpy as np
+import csv
+from netCDF4 import Dataset
+from scipy.interpolate import LinearNDInterpolator
+import matplotlib.pyplot as plt
+from matplotlib.pyplot import cm
+import matplotlib.pyplot as plt
+import matplotlib.tri as tri
+from matplotlib.colorbar import Colorbar
+from matplotlib.colors import Normalize, TwoSlopeNorm
+
+
+def plot_transect(data_path, variable, ax, times=[0],
+                  coords_file=None, x=None, y=None,
+                  time_cmap_name="plasma",
+                  temperature_cmap_name='YlGnBu_r'):
+    """
+    Plot transects of desired variable from MALI output.
+
+    Parameters
+    ----------
+    data_path : str
+        Path to .nc file containing MALI mesh and variables to plot
+
+    variable : str
+        MALI variable to plot. Can also be "geometry", which will
+        calculate upper and lower surfaces from thickness and bed topography.
+
+    ax : matplotlib.axes._axes.Axes
+        Axes on which to plot variable
+
+    times : list of ints, optional
+        Time indices at which to plot variable.
+
+    coords_file : str, optional
+        Path to file containing coordinates for transect
+
+    x : list of floats, optional
+        x coordinates defining transect if not using coords_file
+
+    y : list of floats, optional
+        y coordinates defining transect if not using coords_file
+
+    time_cmap_name : str, optional
+        Name of matplotlib colormap for multiple time levels odataset.variable
+
+    temperature_cmap_nam : str, optional
+        Name of matplotlib colormap for temperature transect
+    """
+
+    dataset = Dataset(data_path)
+    xCell = dataset.variables["xCell"][:]
+    yCell = dataset.variables["yCell"][:]
+    nCells = dataset.dimensions['nCells'].size
+    areaCell = dataset.variables["areaCell"][:]
+    layerThicknessFractions = dataset.variables["layerThicknessFractions"][:]
+    nVertLevels = dataset.dimensions['nVertLevels'].size
+    thk = dataset.variables["thickness"][:]
+    bed = dataset.variables["bedTopography"][:]
+    # ensure that times is a list, otherwise indexing will cause errors
+    times = list(times)
+    if "daysSinceStart" in dataset.variables.keys():
+        use_yrs = True
+        yrs = dataset.variables["daysSinceStart"][:] / 365.
+        times_list = [f'{yrs[i]:.1f}' for i in times]
+    else:
+        use_yrs = False
+        times_list =  [str(i) for i in times]
+
+    # Replace -1 time index with last forward index of time array.
+    # It is unclear why these need to be in separate if-statements, but they do.
+    if -1 in times:
+        times[times.index(-1)] = int(dataset.dimensions['Time'].size - 1)
+    if '-1' in times_list:
+        times_list[times_list.index('-1')] = str(dataset.dimensions['Time'].size - 1)
+
+    # Use coordinates from CSV file or x,y options, but not both.
+    # CSV file takes precedent if both are present.
+    if coords_file is not None:
+        x = []
+        y = []
+        with open(options.coords_file, newline='') as csvfile:
+            reader = csv.reader(csvfile, delimiter=',')
+
+            for row in reader:
+                x.append(float(row[0]))
+                y.append(float(row[1]))
+        if [options.x_coords, options.y_coords] != [None, None]:
+            print('-c and -x/-y options were both provided. Reading from ',
+                  f'{options.coords_file} and ignoring -x and -y settings.')
+        x = np.asarray(x)
+        y = np.asarray(y)
+
+    # increase sampling to match highest mesh resolution
+    total_distance, = np.cumsum( np.sqrt( np.diff(x)**2. + np.diff(y)**2. ) )
+    n_samples = int(round(total_distance / np.min(dataset.variables["dcEdge"][:])))
+    x_interp = np.interp(np.linspace(0, len(x)-1, n_samples),
+                         np.linspace(0, len(x)-1, len(x)), x)
+    y_interp = np.interp(np.linspace(0, len(y)-1, n_samples),
+                         np.linspace(0, len(y)-1, len(y)), y)
+
+    d_distance = np.zeros(len(x_interp))
+    for ii in np.arange(1, len(x_interp)):
+        d_distance[ii] = np.sqrt( (x_interp[ii] - x_interp[ii-1])**2 +
+                                  (y_interp[ii] - y_interp[ii-1])**2 )
+
+    distance = np.cumsum(d_distance) / 1000.  # in km for plotting
+
+    # Handle colormaps
+    time_cmap = plt.get_cmap(time_cmap_name)
+    time_colors = time_cmap(np.linspace(0,1,len(times)))
+
+    if variable in ["temperature", "geometry"]:
+        # Assume constant bed
+        bed_interpolator = LinearNDInterpolator(
+                     np.vstack((xCell, yCell)).T, bed[0,:])
+        bed_transect = bed_interpolator(np.vstack((x_interp, y_interp)).T)
+        ax.plot(distance, bed_transect, color='black')
+        for i, time in enumerate(times):
+            thk_interpolator = LinearNDInterpolator(
+                np.vstack((xCell, yCell)).T, thk[time,:])
+            thk_transect = thk_interpolator(np.vstack((x_interp, y_interp)).T)
+            lower_surf = np.maximum( -910. / 1028. * thk_transect, bed_transect)
+            lower_surf_nan = lower_surf.copy()  # for plotting
+            lower_surf_nan[thk_transect==0.] = np.nan
+            upper_surf = lower_surf + thk_transect
+            upper_surf_nan = upper_surf.copy()  # for plotting
+            upper_surf_nan[thk_transect==0.] = np.nan
+            ax.plot(distance, lower_surf_nan, color=time_colors[i])
+            ax.plot(distance, upper_surf_nan, color=time_colors[i])
+
+        if variable == "temperature":
+            time = times[-1]
+            if len(times) > 1:
+                print("Cannot plot temperature at more than one time."
+                      " Only plotting temperature for final time.")
+            temperature = dataset.variables["temperature"][:]
+            layer_thk = np.zeros((len(thk_transect), nVertLevels + 1))
+            layer_midpoints = np.zeros((len(thk_transect), nVertLevels))
+            layer_interfaces = np.zeros((len(thk_transect), nVertLevels + 1))
+            layer_thk[:,0] = 0.
+            for ii in range(len(thk_transect)):
+                layer_thk[ii,1:] = np.cumsum(layerThicknessFractions *
+                                             thk_transect[ii])
+                layer_midpoints[ii,:] = upper_surf[ii] - (layer_thk[ii,1:] +
+                                                         layer_thk[ii,0:-1]) / 2.
+                layer_interfaces[ii,:] = upper_surf[ii] - layer_thk[ii,:]
+
+            temp_transect = np.zeros((len(x_interp), nVertLevels))
+            for lev in range(nVertLevels):
+                print(f'Interpolating temperature for level {lev} of {nVertLevels}. ', end="")
+                temp_interpolant = LinearNDInterpolator(
+                                        np.vstack((xCell, yCell)).T,
+                                        temperature[time,:,lev])
+                temp_transect[:, lev] = temp_interpolant(
+                                            np.vstack((x_interp, y_interp)).T)
+                temp_transect[temp_transect == 0.] = np.nan
+            temp_plot = ax.pcolormesh(np.tile(distance, (nVertLevels+1,1)).T,
+                                layer_interfaces[:,:], temp_transect[1:,:],
+                                cmap=temperature_cmap_name,
+                                vmin=240., vmax=273.15)
+            temp_cbar = plt.colorbar(temp_plot)
+            temp_cbar.set_label('Temperature (K)')
+            temp_cbar.ax.tick_params(labelsize=12)
+
+    else:
+        var = dataset.variables[variable][:]
+        for i, time in enumerate(times):
+            var_interpolator = LinearNDInterpolator(
+                                  np.vstack((xCell, yCell)).T, var[i,:])
+            var_transect = var_interpolator(np.vstack((x_interp, y_interp)).T)
+            ax.plot(distance, var_transect, color=time_colors[i])
+    if (len(times) > 1):
+        time_cbar = plt.colorbar(cm.ScalarMappable(cmap=time_cmap_name), ax=ax)
+        time_cbar.ax.tick_params(labelsize=12)
+        if use_yrs:
+            time_cbar.set_label('Year')
+        else:
+            time_cbar.set_label('time index')
+        time_cbar.set_ticks(times / np.max(times))
+        time_cbar.set_ticklabels(times_list)
+
+
+def plot_map(data_path, variable, ax, time=0, cmap=None,
+             vmin=None, vmax=None, norm=None, log_plot=False,
+             mesh_file=None, triangles=None,
+             plot_grounding_line=False, variable_name=None):
+    """
+    Plot map of MALI output either by specifying a variable name or
+    a pre-computed field.
+
+    Parameters
+    ----------
+    data_path : str
+        Path to .nc file containing variables to plot. May contain
+        MALI mesh fields, or you can use the mesh_file argument.
+
+    variable : str or numpy.ndarray
+        MALI variable to plot. If a string is specified, the variable with that
+        name will be read from the .nc file at data_path. If a numpy array is
+        given, that array will be plotted directly.
+
+    ax : matplotlib.axes._axes.Axes
+        Axes on which to plot variable
+
+    time : int, optional
+        Time index at which to plot variable.
+
+    cmap : str, optional
+        Name of matplotlib colormap for multiple time levels of variable
+
+    vmin : float, optional
+        Minimum value to use for plotting. If not specified, the 1st
+        percentile value will be used (not weighted by cell area).
+
+    vmax : float, optional
+        Maximum value to use for plotting. If not specified, the 99th
+        percentile value will be used (not weighted by cell area).
+
+    log_plot : boolean, optional
+        Whether to plot log10(variable)
+
+    mesh_file : str, optional
+        Optional file used to specify mesh variables. If not provided, mesh
+        variables will be read from data_path
+
+    triangles : matplotlib.tri._triangulation.Triangulation, optional
+        Triangles to use for plotting. If not specified,
+        they will be calculated.
+
+    plot_grounding_line : boolean, optional
+        Whether to plot the grounding line along with variable.
+
+    variable_name : str
+        Name to use for colorbar if specifying `variable` as a numpy array.
+
+    Returns
+    -------
+    var_plot : matplotlib.collections.PolyCollection
+        Plot of variable
+
+    cbar : matplotlib.colorbar.Colorbar
+        Colorbar object corresponded to var_plot
+
+    gl_plot : matplotlib.tri._tricontour.TriContourSet or None
+        Contour object of the grounding line, if plot_grounding_line=True
+    """
+
+    sec_per_year = 60. * 60. * 24. * 365.
+
+    dataset = Dataset(data_path)
+    dataset.set_auto_mask(False)
+
+    if triangles is None:
+        if mesh_file is None:
+           mesh = Dataset(data_path)
+        else:
+           mesh = Dataset(mesh_file)
+
+        triangles, tri_mask = _get_triangles(mesh)
+        mesh.close()
+
+    if type(variable) is str:
+        variable_name = variable
+        if variable == 'observedSpeed':
+            var_to_plot = np.sqrt(dataset.variables['observedSurfaceVelocityX'][:]**2 +
+                                  dataset.variables['observedSurfaceVelocityY'][:]**2)
+        else:
+            var_to_plot = dataset.variables[variable][:]
+    else:
+        var_to_plot = variable
+
+    if len(np.shape(var_to_plot)) == 1:
+       var_to_plot = var_to_plot.reshape((1, np.shape(var_to_plot)[0]))
+
+    if 'Speed' in variable:
+        units = 'm yr^{-1}'
+        var_to_plot *= sec_per_year
+    else:
+        try:
+            units = f'({dataset.variables[variable].units})'
+        except AttributeError:
+            units = "{}"  # This leaves out units on the colorbar
+        except TypeError:
+            units = "{}"
+
+    default_colors = {'thickness' : 'Blues',
+                     'surfaceSpeed' : 'plasma',
+                     'basalSpeed' : 'plasma',
+                     'bedTopography' : 'BrBG',
+                     'floatingBasalMassBalApplied' : 'cividis'
+                    }
+    # List of diverging colormaps for use in plotting bedTopography.
+    # I don't see a way around hard-coding this.
+    div_color_maps = ['PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu', 'RdYlBu',
+                          'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic']
+
+    if cmap is None:
+        if variable_name in default_colors.keys():
+            cmap = default_colors[variable]
+        else:
+            cmap = "viridis"
+
+    if log_plot:
+        var_to_plot = np.log10(var_to_plot)
+        # Get rid of +/- inf values that ruin vmin and vmax
+        # calculations below.
+        var_to_plot[np.isinf(var_to_plot)] = np.nan
+        colorbar_label_prefix = 'log10 '
+    else:
+        colorbar_label_prefix = ''
+
+    # Set lower and upper bounds for plotting
+    if vmin is None:
+        # 0.1 m/yr is a pretty good lower bound for speed
+        first_quant = np.nanquantile(var_to_plot[time, :], 0.01)
+        if 'Speed' in variable and log_plot:
+            vmin = max(first_quant, -1.)
+        else:
+            vmin = first_quant
+    if vmax is None:
+        vmax = np.nanquantile(var_to_plot[time, :], 0.99)
+    # Plot bedTopography on an asymmetric colorbar if appropriate
+    if norm is None:
+        if ( (variable_name == 'bedTopography') and
+             (np.nanquantile(var_to_plot[time, :], 0.99) > 0.) and
+             (cmap in div_color_maps) ):
+            norm = TwoSlopeNorm(vmin=vmin, vmax=vmax, vcenter=0.)
+        else:
+            norm = Normalize(vmin=vmin, vmax=vmax)
+
+    var_plot = ax.tripcolor(
+        triangles, var_to_plot[time, :], cmap=cmap,
+        shading='flat', norm=norm)
+    ax.set_aspect('equal')
+
+    cbar = plt.colorbar(ax=ax, mappable=var_plot,
+                        orientation='vertical',
+                        label=f'{colorbar_label_prefix}{variable_name} ${units}$')
+    if plot_grounding_line:
+        valid_masks, grounding_line_mask, _, _, _ = _calculate_masks(dataset)
+        if valid_masks:
+            gl_plot = ax.tricontour(triangles, grounding_line_mask[time, :],
+                                    levels=[0.9999], colors='white',
+                                    linestyles='solid')
+        else:
+            gl_plot = None
+    else:
+        gl_plot = None
+
+    return var_plot, cbar, gl_plot
+
+
+def plot_grounding_lines(data_paths, ax, times=[0],
+                         cmap="plasma_r", norm=None,
+                         mesh_file=None, triangles=None):
+    """
+    Plot MALI grounding line at arbitrary number of times.
+
+    Parameters
+    ----------
+    data_paths : str or list of str
+        Path(s) to MALI file. May contain MALI mesh fields,
+        or you can use the mesh_file argument.
+
+    ax : matplotlib.axes._axes.Axes
+        Axes on which to plot variable
+
+    time : list of ints, optional
+        Time indices at which to plot variable.
+
+    cmap : str, optional
+        Name of matplotlib colormap for multiple time levels of variable
+
+    mesh_file : str, optional
+        Optional file used to specify mesh variables. If not provided, mesh
+        variables will be read from data_path
+
+    triangles : matplotlib.tri._triangulation.Triangulation, optional
+        Triangles to use for plotting. If not specified,
+        they will be calculated.
+
+    Returns
+    -------
+    gl_plots : list of matplotlib.tri._tricontour.TriContourSet
+        List of grounding line contour objects
+    """
+
+    # Ensure data_path is a list for flexibility in loops
+    if type(data_paths) != list:
+        data_paths = [data_paths]
+
+    # If triangles are not specified, use first
+    # data file to define mesh, or use mesh_file
+    if triangles is None:
+        if mesh_file is None:
+            mesh = Dataset(data_paths[0])
+        else:
+            mesh = Dataset(mesh_file)
+
+        triangles, tri_mask = _get_triangles(mesh)
+        mesh.close()
+
+    # Loop over all files and time levels to create
+    # lists of the grounding lines, and their associated years
+    plot_times = []
+    grounding_line_masks = []
+    gl_plots = []
+    for file in data_paths:
+        f = Dataset(file)
+        f.set_auto_mask(False)
+        if 'daysSinceStart' in f.variables.keys():
+            yr = f.variables['daysSinceStart'][times] / 365.
+        else:
+            yr = times
+        plot_times.append(yr)
+        valid_masks, grounding_line_mask, _, _, _ = _calculate_masks(f)
+        if valid_masks:
+            for time in times:
+                grounding_line_masks.append(grounding_line_mask[time, :])
+        f.close()
+
+    # Determine mapping between plot time and colormap.
+    # If just plotting one time, then use maximum value
+    # of the specified colormap.
+    n_times = len(times) * len(data_paths)
+    gl_cmap = plt.get_cmap(cmap)
+    if n_times > 1:
+        plot_times = np.squeeze(np.ravel(plot_times))
+        plot_times_norm = ( (plot_times - np.min(plot_times)) /
+                           np.max(plot_times - np.min(plot_times)) )
+    else:
+        plot_times_norm = np.ones_like(plot_times)
+
+    time_colors = gl_cmap(plot_times_norm)
+
+    for ii, mask in enumerate(grounding_line_masks):
+        gl_plots.append(ax.tricontour(triangles, mask,
+                                levels=[0.9999], linestyles='solid',
+                                colors=time_colors[ii, None]))
+
+    if len(plot_times) > 1:
+        time_cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax,
+                                 location='bottom', label="Grounding line year")
+        time_cbar.ax.tick_params(labelsize=12)
+        time_cbar.set_ticks(plot_times_norm)
+        time_cbar.set_ticklabels(str(i) for i in plot_times)
+    else:
+        time_cbar = None
+
+    return gl_plots, time_cbar
+
+
+def _dist(i1, i2, xCell, yCell):
+
+    dist = ((xCell[i1]-xCell[i2])**2 + (yCell[i1]-yCell[i2])**2)**0.5
+    return dist
+
+
+def _calculate_masks(dataset):
+
+    # Set bitmask values
+    initial_extent_value = 1
+    dynamic_value = 2
+    float_value = 4
+    grounding_line_value = 256
+    rhoi = 910.
+    rhosw = 1028.
+
+    if 'cellMask' in dataset.variables.keys():
+        valid_masks = True
+        cellMask = dataset.variables["cellMask"][:]
+        float_mask = (cellMask & float_value) // float_value
+        dynamic_mask = (cellMask & dynamic_value) // dynamic_value
+        grounding_line_mask = (cellMask & grounding_line_value) // grounding_line_value
+        initial_extent_mask = (cellMask & initial_extent_value) // initial_extent_value
+    elif ( 'cellMask' not in dataset.variables.keys() and
+         'thickness' in dataset.variables.keys() and
+         'bedTopography' in dataset.variables.keys() ):
+        print(f'cellMask is not present in output file {run};'
+               ' calculating masks from ice thickness')
+        valid_masks = True
+        grounded_mask = (dataset.variables['thickness'][:] >
+                        (-rhosw / rhoi *
+                         dataset.variables['bedTopography'][:]))
+        # This isn't technically correct, but works for plotting
+        grounding_line_mask = grounded_mask.copy()
+        initial_extent_mask = (dataset.variables['thickness'][:] > 0.)
+    else:
+        print('cellMask and thickness and/or bedTopography'
+              f' not present in output file {run};'
+               ' Skipping mask calculation.')
+        valid_masks = False
+
+    return valid_masks, grounding_line_mask, float_mask, dynamic_mask, initial_extent_mask
+
+
+def _get_triangles(mesh):
+
+    xCell = mesh.variables["xCell"][:]
+    yCell = mesh.variables["yCell"][:]
+    dcEdge = mesh.variables["dcEdge"][:]
+
+    triang = tri.Triangulation(xCell, yCell)
+    tri_mask = np.zeros(len(triang.triangles))
+
+    # Maximum distance in m of edges between points.
+    # Make twice dcEdge to be safe
+    max_dist = np.max(dcEdge) * 2.0
+    for t in range(len(triang.triangles)):
+        thisTri = triang.triangles[t, :]
+        if _dist(thisTri[0], thisTri[1], xCell, yCell) > max_dist:
+            tri_mask[t] = True
+        if _dist(thisTri[1], thisTri[2], xCell, yCell) > max_dist:
+            tri_mask[t] = True
+        if _dist(thisTri[0], thisTri[2], xCell, yCell) > max_dist:
+            tri_mask[t] = True
+    triang.set_mask(tri_mask)
+
+    return triang, tri_mask