diff --git a/config.analysis b/config.analysis
index 4605721ff..99a1b489c 100644
--- a/config.analysis
+++ b/config.analysis
@@ -82,6 +82,9 @@ generate = 1
 [seaice_modelvsobs]
 generate = 1
 
+[moc_postprocess]
+generate = 0
+
 [ohc_timeseries]
 ## compare to output from another model run?
 #compare_with_model = True
@@ -131,6 +134,18 @@ compare_with_obs = True
 # output, N_movavg=12 corresponds to a 12-month moving average window)
 N_movavg = 12
 
+
+[moc_postprocess]
+filename  = ../mpaso.hist.am.timeSeriesStats.yr100-119.avg.nc
+filenamemesh = ../mpaso.rst.0002-01-01_00000.nc
+filenamemask = ../EC60to30v3_SingleRegionAtlanticWTransportTransects_masks.nc
+timeIndex = 0
+horVelName = time_avg_normalVelocity
+vertVelName = time_avg_vertVelocityTop
+contourLines = np.arange(-8,20.1,2)
+# This is specific to Atlantic MOC:
+mocLat = np.arange(-34.5,70.1,0.5)
+
 [seaice_timeseries]
 ## compare to output from another model run?
 #compare_with_model = True
diff --git a/mpas_analysis/configuration/MpasAnalysisConfigParser.py b/mpas_analysis/configuration/MpasAnalysisConfigParser.py
index fe1cd7ddf..8a73451d0 100644
--- a/mpas_analysis/configuration/MpasAnalysisConfigParser.py
+++ b/mpas_analysis/configuration/MpasAnalysisConfigParser.py
@@ -1,16 +1,17 @@
-import numbers
-import ast
-
 """
 A configuratin parser class for MPAS analysis.  MpasAnalysisConfigParser adds
 the capabilities to get an option including a default value
 (`getWithDefault(section, option, default, ...)`) and to get options
 that are lists, tuples, dicts, etc (`getExpression(section, option)`).
 
-Author: Xylar Asay-Davis
-Last Modified: 12/07/2016
+Author: Xylar Asay-Davis, Phillip J. Wolfram
+Last Modified: 01/31/2017
 """
 
+import numbers
+import ast
+import numpy as np
+import numpy
 from ConfigParser import ConfigParser
 
 
@@ -44,7 +45,7 @@ def getWithDefault(self, section, option, default):
         self.set(section, option, str(default))
         return default
 
-    def getExpression(self, section, option, elementType=None):
+    def getExpression(self, section, option, elementType=None, usenumpy=False):
         """
         Get an option as an expression (typically a list, though tuples and
         dicts should also work).  `section` and `option` work as in `get(...)`.
@@ -56,11 +57,17 @@ def getExpression(self, section, option, elementType=None):
         useful for ensuring that all elements in a list of numbers are of type
         float, rather than int, when the distinction is important.
 
-        Author: Xylar Asay-Davis
-        Last Modified: 12/0y/2016
+        If `numpy` is True, expression is evaluated within the context of
+        having numpy and / or np available.
+
+        Author: Xylar Asay-Davis, Phillip J. Wolfram
+        Last Modified: 01/31/2017
         """
         expressionString = self.get(section, option)
-        result =  ast.literal_eval(expressionString)
+        if usenumpy:
+            result = eval(expressionString)
+        else:
+            result =  ast.literal_eval(expressionString)
 
         if elementType is not None:
             if isinstance(result, (list, tuple)):
diff --git a/mpas_analysis/ocean/meridional_overturning_circulation.py b/mpas_analysis/ocean/meridional_overturning_circulation.py
new file mode 100644
index 000000000..a6497685c
--- /dev/null
+++ b/mpas_analysis/ocean/meridional_overturning_circulation.py
@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+"""
+Compute and plot moc, as a post-processing computation
+
+Author: Mark Petersen, Phillip J. Wolfram
+Last Modified: 01/20/2017
+"""
+
+from netCDF4 import Dataset
+import numpy as np
+from ..shared.constants.constants import m3ps_to_Sv
+from ..shared.plot.plotting import plot_moc
+
+def moc_postprocess(config, streamMap=None, variableMap=None): #{{{
+    """
+    Plots a comparison of ACME/MPAS output to SST or MLD observations
+
+    config is an instance of MpasAnalysisConfigParser containing configuration
+    options.
+
+    If present, streamMap is a dictionary of MPAS-O stream names that map to
+    their mpas_analysis counterparts.
+
+    If present, variableMap is a dictionary of MPAS-O variable names that map
+    to their mpas_analysis counterparts.
+
+    Authors: Mark Petersen, Phillip J. Wolfram
+    Modified: 01/20/2017
+    """
+
+    # read parameters from config file
+    fileName = config.get('moc_postprocess', 'filename')
+    fileNameMesh = config.get('moc_postprocess', 'filenamemesh')
+    fileNameMask = config.get('moc_postprocess', 'filenamemask')
+    timeIndex = config.get('moc_postprocess', 'timeIndex')
+    horVelName = config.get('moc_postprocess', 'horVelName')
+    vertVelName = config.get('moc_postprocess', 'vertVelname')
+    mocLat = config.getExpression('moc_postprocess', 'mocLat', usenumpy=True)
+
+    # open files, load variables
+    print('** load variables: \n')
+
+    print 'Reading data from: ' + fileName
+    ncfile = Dataset(fileName,'r')
+    nCells = ncfile.dimensions['nCells'].size
+    nVertLevels = ncfile.dimensions['nVertLevels'].size
+
+    horVel = ncfile.variables[horVelName][:,:,:]
+    vertVel = ncfile.variables[vertVelName][:,:,:]
+
+    print 'Reading data from: ' + fileNameMesh
+    ncfileMesh = Dataset(fileNameMesh,'r')
+    dvEdge = ncfileMesh.variables['dvEdge'][:]
+    areaCell = ncfileMesh.variables['areaCell'][:]
+    latCell = ncfileMesh.variables['latCell'][:]*180./np.pi
+    refBottomDepth = ncfileMesh.variables['refBottomDepth'][:]
+
+    refLayerThickness = np.zeros(nVertLevels)
+    refLayerThickness[0] = refBottomDepth[0]
+    refLayerThickness[1:nVertLevels] = refBottomDepth[1:nVertLevels]-refBottomDepth[0:nVertLevels-1]
+
+    refTopDepth = np.zeros(nVertLevels+1)
+    refTopDepth[1:nVertLevels+1] = refBottomDepth[0:nVertLevels]
+
+    print 'Reading transects from: ' + fileNameMask
+    ncfileMask = Dataset(fileNameMask,'r')
+    nTransects = ncfileMask.dimensions['nTransects'].size
+    maxEdgesInTransect = ncfileMask.dimensions['maxEdgesInTransect'].size
+
+    transectEdgeMaskSigns = ncfileMask.variables['transectEdgeMaskSigns'][:,:]
+    transectEdgeGlobalIDs = ncfileMask.variables['transectEdgeGlobalIDs'][:,:]
+
+    print 'Reading regionss from: ' + fileNameMask
+    ncfileMask = Dataset(fileNameMask,'r')
+    nRegions = ncfileMask.dimensions['nRegions'].size
+    regionCellMasks = ncfileMask.variables['regionCellMasks']
+
+    # compute transport through transects
+    print('** compute transport: \n')
+
+    # the volume transport
+    transport = np.zeros(nTransects);
+    transportZ = np.zeros([nVertLevels,nTransects]);
+    transportZEdge = np.zeros([nVertLevels,maxEdgesInTransect,nTransects]);
+
+    for iTransect in range(nTransects):
+        for i in range(maxEdgesInTransect):
+            if transectEdgeGlobalIDs[iTransect, i]==0:
+                break
+            iEdge = transectEdgeGlobalIDs[iTransect, i] - 1 # subtract 1 because of python 0-indexing
+            for k in range(nVertLevels):
+                transportZEdge[k,i,iTransect] = ( transectEdgeMaskSigns[iEdge,iTransect]
+                        * horVel[timeIndex,iEdge,k] * dvEdge[iEdge] * refLayerThickness[k]*m3ps_to_Sv )
+                transportZ[k,iTransect] = transportZ[k,iTransect] + transportZEdge[k,i,iTransect]
+                transport[iTransect] = transport[iTransect] + transportZEdge[k,i,iTransect]
+
+    # compute MOC
+    print('** compute moc: \n')
+
+    nLat = np.size(mocLat)
+    mocTop = np.zeros((nLat,nVertLevels+1))
+
+    for iRegion in range(nRegions):
+        # assume transects and regions have the same index ordering:
+       iTransect = iRegion
+
+       for k in range(1,nVertLevels+1):
+           mocTop[0,k] = mocTop[0,k-1] + transportZ[k-1,iTransect]/m3ps_to_Sv
+
+       for iLat in range(1,nLat):
+           ind =np.logical_and(np.logical_and(regionCellMasks[:,iRegion]==1,latCell >= mocLat[iLat-1] ), latCell <  mocLat[iLat])
+
+           for k in range(1,nVertLevels+1):
+               mocTop[iLat,k] = mocTop[iLat-1,k] + np.sum(vertVel[timeIndex,ind,k]*areaCell[ind])
+
+       # convert m^3/s to Sverdrup
+       mocTop = mocTop * m3ps_to_Sv
+
+    # plot MOC
+    print('** plot moc: \n')
+    plot_moc(config, mocLat, refTopDepth, mocTop)
+
+    ncfile.close()
+    ncfileMask.close()
+
+    return #}}}
+
+# vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python
diff --git a/mpas_analysis/shared/constants/constants.py b/mpas_analysis/shared/constants/constants.py
index 1e1a72f6d..c08fe2cfd 100644
--- a/mpas_analysis/shared/constants/constants.py
+++ b/mpas_analysis/shared/constants/constants.py
@@ -22,4 +22,5 @@
 
 dinmonth = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
 
+m3ps_to_Sv = 1e-6; # m^3/sec flux to Sverdrups
 
diff --git a/mpas_analysis/shared/plot/plotting.py b/mpas_analysis/shared/plot/plotting.py
index ee72a1834..9fd232eb1 100644
--- a/mpas_analysis/shared/plot/plotting.py
+++ b/mpas_analysis/shared/plot/plotting.py
@@ -1,12 +1,27 @@
+#!/usr/bin/env python
+"""
+Plotting routines
+
+Xylar Asay-Davis, Phillip Wolfram
+Last modified: 01/23/2017
+"""
+
 import matplotlib.pyplot as plt
 import pandas as pd
 from mpl_toolkits.basemap import Basemap
 import matplotlib.colors as cols
 import numpy as np
 
-def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, fileout,
-                             lineStyles, lineWidths, title_font_size=None,
-                             figsize=(15,6), dpi=300):
+def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel,
+                             fileout,lineStyles, lineWidths,
+                             title_font_size=None,
+                             figsize=(15, 6), dpi=300): # {{{
+    """
+    Timeseries plotting routine.
+
+    Xylar Asay-Davis
+    Last modified: 09/15/2016
+    """
 
     plt.figure(figsize=figsize, dpi=dpi)
 
@@ -31,11 +46,13 @@ def timeseries_analysis_plot(config, dsvalues, N, title, xlabel, ylabel, fileout
     if (ylabel != None):
         plt.ylabel(ylabel, **axis_font)
     if (fileout is not None):
-        plt.savefig(fileout,dpi=dpi,bbox_inches='tight',pad_inches=0.1)
+        plt.savefig(fileout,dpi=dpi, bbox_inches='tight', pad_inches=0.1)
 
-    if not config.getboolean('plot','displayToScreen'):
+    if not config.getboolean('plot', 'displayToScreen'):
       plt.close()
 
+    return # }}}
+
 
 def plot_polar_comparison(
     config,
@@ -58,8 +75,14 @@ def plot_polar_comparison(
     diffTitle = "Model-Observations",
     cbarlabel = "units",
     title_font_size = None,
-    figsize = (8,22),
-    dpi = 300):
+    figsize = (8, 22),
+    dpi = 300): # {{{
+    """
+    Plots model versus observations for polar regions.
+
+    Xylar Asay-Davis
+    Last modified: 09/15/2016
+    """
 
     # set up figure
     fig = plt.figure(figsize=figsize, dpi=dpi)
@@ -72,54 +95,61 @@ def plot_polar_comparison(
         fig.suptitle(title, y=0.95, **title_font)
     axis_font = {'size':config.get('plot', 'axis_font_size')}
 
-    m = Basemap(projection=plotProjection,boundinglat=latmin,lon_0=lon0,resolution='l')
+    m = Basemap(projection=plotProjection, boundinglat=latmin,
+                lon_0=lon0, resolution='l')
     x, y = m(Lons, Lats) # compute map proj coordinates
 
     normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N)
     normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N)
 
-    plt.subplot(3,1,1)
+    plt.subplot(3, 1, 1)
     plt.title(modelTitle, y=1.06, **axis_font)
     m.drawcoastlines()
     m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,81.,10.))
-    m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True])
-    cs = m.contourf(x,y,modelArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs)
-    cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsModelObs, boundaries=clevsModelObs)
+    m.drawparallels(np.arange(-80., 81., 10.))
+    m.drawmeridians(np.arange(-180., 181., 20.),
+                    labels=[True, True, True, True])
+    cs = m.contourf(x, y, modelArray, cmap=cmapModelObs, norm=normModelObs,
+                    spacing='uniform', levels=clevsModelObs)
+    cbar = m.colorbar(cs, location='right', pad="15%", spacing='uniform',
+                      ticks=clevsModelObs, boundaries=clevsModelObs)
+
     cbar.set_label(cbarlabel)
 
-    plt.subplot(3,1,2)
+    plt.subplot(3, 1, 2)
     plt.title(obsTitle, y=1.06, **axis_font)
     m.drawcoastlines()
     m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,81.,10.))
-    m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True])
-    cs = m.contourf(x,y,obsArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs)
-    cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsModelObs, boundaries=clevsModelObs)
+    m.drawparallels(np.arange(-80., 81., 10.))
+    m.drawmeridians(np.arange(-180., 181., 20.),
+                    labels=[True, True, True, True])
+    cs = m.contourf(x, y, obsArray, cmap=cmapModelObs, norm=normModelObs,
+                    spacing='uniform', levels=clevsModelObs)
+    cbar = m.colorbar(cs, location='right', pad="15%", spacing='uniform',
+                      ticks=clevsModelObs, boundaries=clevsModelObs)
     cbar.set_label(cbarlabel)
 
-    plt.subplot(3,1,3)
+    plt.subplot(3, 1, 3)
     plt.title(diffTitle, y=1.06, **axis_font)
     m.drawcoastlines()
     m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,81.,10.))
-    m.drawmeridians(np.arange(-180.,181.,20.),labels=[True,True,True,True])
-    cs = m.contourf(x,y,diffArray,cmap=cmapDiff,norm=normDiff,spacing='uniform',levels=clevsDiff)
-    cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',ticks=clevsDiff,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="15%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsDiff, boundaries=clevsDiff)
+    m.drawparallels(np.arange(-80., 81., 10.))
+    m.drawmeridians(np.arange(-180., 181., 20.),
+                    labels=[True, True, True, True])
+    cs = m.contourf(x, y, diffArray, cmap=cmapDiff, norm=normDiff,
+                    spacing='uniform', levels=clevsDiff)
+    cbar = m.colorbar(cs,location='right', pad="15%", spacing='uniform',
+                      ticks=clevsDiff, boundaries=clevsModelObs)
     cbar.set_label(cbarlabel)
 
     if (fileout is not None):
-        plt.savefig(fileout,dpi=dpi,bbox_inches='tight',pad_inches=0.1)
+        plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1)
 
-    if not config.getboolean('plot','displayToScreen'):
+    if not config.getboolean('plot', 'displayToScreen'):
         plt.close()
 
+    return # }}}
+
 def plot_global_comparison(
     config,
     Lons,
@@ -138,8 +168,14 @@ def plot_global_comparison(
     diffTitle = "Model-Observations",
     cbarlabel = "units",
     title_font_size = None,
-    figsize = (8,12),
-    dpi = 300):
+    figsize = (8, 12),
+    dpi = 300): # {{{
+    """
+    Plots model versus observations for global projections.
+
+    Xylar Asay-Davis
+    Last modified: 09/15/2016
+    """
 
     # set up figure
     fig = plt.figure(figsize=figsize, dpi=dpi)
@@ -151,54 +187,102 @@ def plot_global_comparison(
                       'weight':config.get('plot', 'title_font_weight')}
         fig.suptitle(title, y=0.95, **title_font)
     axis_font = {'size':config.get('plot', 'axis_font_size')}
-        
-    m = Basemap(projection='cyl',llcrnrlat=-85,urcrnrlat=86,llcrnrlon=-180,urcrnrlon=181,resolution='l')
-    #m = Basemap(projection='robin',lon_0=200,resolution='l') # this doesn't work because lons are -180 to 180..
+
+    m = Basemap(projection='cyl', llcrnrlat=-85, urcrnrlat=86, llcrnrlon=-180,
+                urcrnrlon=181, resolution='l')
+    # this doesn't work because lons are -180 to 180..
+    #m = Basemap(projection='robin', lon_0=200, resolution='l')
     x, y = m(Lons, Lats) # compute map proj coordinates
 
     normModelObs = cols.BoundaryNorm(clevsModelObs, cmapModelObs.N)
     normDiff = cols.BoundaryNorm(clevsDiff, cmapDiff.N)
-        
-    plt.subplot(3,1,1)
+
+    plt.subplot(3, 1, 1)
     plt.title(modelTitle, y=1.06, **axis_font)
     m.drawcoastlines()
-    m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,80.,20.),labels=[True,False,False,False])
-    m.drawmeridians(np.arange(-180.,180.,60.),labels=[False,False,False,True])
-    cs = m.contourf(x,y,modelArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs,extend='both')
-    cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsModelObs, boundaries=clevsModelObs)
+    m.fillcontinents(color='grey', lake_color='white')
+    m.drawparallels(np.arange(-80., 80., 20.),
+                    labels=[True, False, False, False])
+    m.drawmeridians(np.arange(-180., 180., 60.),
+                    labels=[False, False, False, True])
+    cs = m.contourf(x, y, modelArray, cmap=cmapModelObs, norm=normModelObs,
+                    spacing='uniform', levels=clevsModelObs, extend='both')
+    cbar = m.colorbar(cs,location='right', pad="5%", spacing='uniform',
+                      ticks=clevsModelObs, boundaries=clevsModelObs)
     cbar.set_label(cbarlabel)
 
-    plt.subplot(3,1,2)
+    plt.subplot(3, 1, 2)
     plt.title(obsTitle, y=1.06, **axis_font)
     m.drawcoastlines()
     m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,80.,20.),labels=[True,False,False,False])
-    m.drawmeridians(np.arange(-180.,180.,40.),labels=[False,False,False,True])   
-    cs = m.contourf(x,y,obsArray,cmap=cmapModelObs,norm=normModelObs,spacing='uniform',levels=clevsModelObs,extend='both')
-    cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsModelObs,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsModelObs, boundaries=clevsModelObs)
+    m.drawparallels(np.arange(-80., 80., 20.),
+                    labels=[True, False, False, False])
+    m.drawmeridians(np.arange(-180., 180., 40.),
+                    labels=[False, False, False, True])
+    cs = m.contourf(x, y, obsArray, cmap=cmapModelObs, norm=normModelObs,
+                    spacing='uniform', levels=clevsModelObs, extend='both')
+    cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform',
+                      ticks=clevsModelObs, boundaries=clevsModelObs)
     cbar.set_label(cbarlabel)
-    
-    plt.subplot(3,1,3)
+
+    plt.subplot(3, 1, 3)
     plt.title(diffTitle, y=1.06, **axis_font)
     m.drawcoastlines()
     m.fillcontinents(color='grey',lake_color='white')
-    m.drawparallels(np.arange(-80.,80.,20.),labels=[True,False,False,False])
-    m.drawmeridians(np.arange(-180.,180.,40.),labels=[False,False,False,True])    
-    cs = m.contourf(x,y,diffArray,cmap=cmapDiff,norm=normDiff,spacing='uniform',levels=clevsDiff,extend='both')
-    cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',ticks=clevsDiff,boundaries=clevsModelObs)
-    #cbar = m.colorbar(cs,location='right',pad="5%",spacing='uniform',extendfrac='auto',
-    #                  extendrect='True',ticks=clevsDiff, boundaries=clevsDiff)
+    m.drawparallels(np.arange(-80., 80., 20.),
+                    labels=[True, False, False, False])
+    m.drawmeridians(np.arange(-180., 180., 40.),
+                    labels=[False, False, False, True])
+    cs = m.contourf(x, y, diffArray, cmap=cmapDiff, norm=normDiff,
+                    spacing='uniform', levels=clevsDiff, extend='both')
+    cbar = m.colorbar(cs, location='right', pad="5%", spacing='uniform',
+                      ticks=clevsDiff, boundaries=clevsModelObs)
     cs.cmap.set_over((1., 1., 1.))
     cs.cmap.set_under((0., 0., 0.))
     cbar.set_label(cbarlabel)
 
     if (fileout is not None):
-        plt.savefig(fileout,dpi=dpi,bbox_inches='tight',pad_inches=0.1)
+        plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1)
 
-    if not config.getboolean('plot','displayToScreen'):
+    if not config.getboolean('plot', 'displayToScreen'):
         plt.close()
+
+    return # }}}
+
+def plot_moc(config, mocLat,
+             refTopDepth,
+             mocTop,
+             fileout = 'moc' + '.png',
+             figsize=(8, 4.5),
+             dpi=300): #{{{
+    """
+    Plots comparison of moc section
+
+    Mark Petersen, Phillip Wolfram
+    Last modified: 01/23/2017
+    """
+
+    contourLines = config.getExpression('moc_postprocess',
+                                        'contourLines', usenumpy=True)
+
+    fig = plt.figure(figsize=figsize, dpi=dpi)
+    X, Y = np.meshgrid(mocLat, refTopDepth) # change to zMid
+
+    contourSet = plt.contourf(X, Y, mocTop.T, levels=contourLines)
+    contourSet2 = plt.contour(X, Y, mocTop.T, levels=contourLines[::2],
+                              colors='k', hold='on')
+
+    plt.xlabel('latitude')
+    plt.ylabel('depth, m')
+    plt.title('MPAS-Ocean Atlantic MOC, Sv')
+    plt.gca().invert_yaxis()
+    cbar = plt.colorbar(contourSet)
+    cbar.add_lines(contourSet2)
+
+    if (fileout is not None):
+        plt.savefig(fileout, dpi=dpi, bbox_inches='tight', pad_inches=0.1)
+
+    if not config.getboolean('plot', 'displayToScreen'):
+        plt.close()
+
+    return #}}}
diff --git a/mpas_analysis/test/test_mpas_config_parser.py b/mpas_analysis/test/test_mpas_config_parser.py
index 4cb09cb1c..b416ecc80 100644
--- a/mpas_analysis/test/test_mpas_config_parser.py
+++ b/mpas_analysis/test/test_mpas_config_parser.py
@@ -1,14 +1,15 @@
 """
 Unit test infrastructure for MpasAnalysisConfigParser
 
-Xylar Asay-Davis
-12/05/2016
+Xylar Asay-Davis, Phillip J. Wolfram
+01/31/2017
 """
 
 import pytest
 from mpas_analysis.test import TestCase, loaddatadir
 from mpas_analysis.configuration.MpasAnalysisConfigParser \
     import MpasAnalysisConfigParser
+from . import requires_numpy
 
 
 @pytest.mark.usefixtures("loaddatadir")
@@ -64,5 +65,15 @@ def test_read_config(self):
                                     'key2': -12,
                                     'key3': False})
 
+    @requires_numpy
+    def test_read_config_numpy(self):
+        self.setup_config()
+
+        # tests numpy evaluation capability
+        import numpy as np
+        result = np.arange(0, 1, 10)
+        for testNumpy in ['testNumpy1', 'testNumpy2']:
+            self.assertArrayEqual(self.config.getExpression('TestNumpy',
+                                              testNumpy, usenumpy=True), result)
 
 # vim: foldmethod=marker ai ts=4 sts=4 et sw=4 ft=python
diff --git a/mpas_analysis/test/test_mpas_config_parser/config.analysis b/mpas_analysis/test/test_mpas_config_parser/config.analysis
index 18071dae5..43ce21f67 100644
--- a/mpas_analysis/test/test_mpas_config_parser/config.analysis
+++ b/mpas_analysis/test/test_mpas_config_parser/config.analysis
@@ -32,3 +32,6 @@ testDict = {'key1': 'string',
             'key2': -12,
             'key3': False}
 
+[TestNumpy]
+testNumpy1 = numpy.arange(0, 1, 10)
+testNumpy2 = np.arange(0, 1, 10)
diff --git a/run_analysis.py b/run_analysis.py
index c02848818..2599a1bb6 100755
--- a/run_analysis.py
+++ b/run_analysis.py
@@ -148,6 +148,12 @@ def analysis(config):  # {{{
         ocn_modelvsobs(config, 'sss', streamMap=oceanStreamMap,
                        variableMap=oceanVariableMap)
 
+    if config.getboolean('moc_postprocess', 'generate'):
+        print ""
+        print "Plotting post-process computed meridional overturning circulation..."
+        from mpas_analysis.ocean.meridional_overturning_circulation import moc_postprocess
+        moc_postprocess(config, streamMap=oceanStreamMap,
+                          variableMap=oceanVariableMap)
 
     # GENERATE SEA-ICE DIAGNOSTICS
     if config.getboolean('seaice_timeseries', 'generate'):