diff --git a/config.analysis b/config.analysis
index 0454e1df2..f664ff72a 100644
--- a/config.analysis
+++ b/config.analysis
@@ -21,6 +21,7 @@ plots_dir = /global/project/projectdirs/acme/xylar/coupled_diagnostics_20160805v
 log_dir = /global/project/projectdirs/acme/xylar/coupled_diagnostics_20160805v0atm.A_WCYCL1850_v0atm.ne30_oEC.edison.alpha7_00-20160520.A_WCYCL1850.ne30_oEC.edison.alpha6_01.logs
 obs_ocndir = /global/project/projectdirs/acme/observations/Ocean
 obs_sstdir = /global/project/projectdirs/acme/observations/Ocean/SST
+obs_sssdir = /global/project/projectdirs/acme/observations/Ocean/SSS
 obs_mlddir = /global/project/projectdirs/acme/observations/Ocean/MLD
 obs_seaicedir = /global/project/projectdirs/acme/observations/SeaIce
 ref_archive_v0_ocndir = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ocn/postprocessing
@@ -70,6 +71,9 @@ generate = 0
 [sst_modelvsobs]
 generate = 1
 
+[sss_modelvsobs]
+generate = 1
+
 [mld_modelvsobs]
 generate = 1
 
@@ -161,6 +165,28 @@ clevsDiff = [-5, -3, -2, -1, 0, 1, 2, 3, 5]
 # Times for comparison times (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, JFM, AMJ, JAS, OND, ANN)
 comparisonTimes =  ['JFM', 'JAS', 'ANN']
 
+[sss_modelvsobs]
+# colormap for model/observations
+#cmapModelObs = viridis
+cmapModelObs = RdYlBu_r
+# colormap for differences
+#cmapDiff = RdBu_r
+cmapDiff = coolwarm
+
+# indices into cmapModelObs for contour color
+cmapIndicesModelObs = [0,40,80,110,140,170,200,230,255]
+# indices into cmapModelObs for contour color
+cmapIndicesDiff = [0,40,80,120,140,170,210,255]
+#cmapIndicesDiff = [0,40,80,127,170,210,255] # good for RdBu_r
+
+# colormap levels/values for contour boundaries
+clevsModelObs = [28,29,30,31,32,33,34,35,36,38]
+clevsDiff = [-3,-2,-1,-0.5,0.5,1,2,3]
+#clevsDiff = [-3,-2,-1,-0.5,0.5,1,2,3] # good for RdBu_r
+
+# Times for comparison times (Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec,JFM,AMJ,JAS,OND,ANN)
+comparisonTimes =  ['JFM','JAS','ANN']
+
 [mld_modelvsobs]
 # colormap for model/observations
 cmapModelObs = viridis
diff --git a/mpas_analysis/ocean/ocean_modelvsobs.py b/mpas_analysis/ocean/ocean_modelvsobs.py
index c3a50cb1b..a459c7175 100644
--- a/mpas_analysis/ocean/ocean_modelvsobs.py
+++ b/mpas_analysis/ocean/ocean_modelvsobs.py
@@ -133,6 +133,35 @@ def ocn_modelvsobs(config, field, streamMap=None, variableMap=None):
         fileOutLabel = "sstHADOI"
         unitsLabel = r'$^o$C'
 
+    elif field == 'sss':
+
+        selvals = {'nVertLevels': 0}
+
+        obs_filename = "{}/Aquarius_V3_SSS_Monthly.nc".format(obsdir)
+        dsData = xr.open_mfdataset(obs_filename)
+
+        time_start = datetime.datetime(2011, 8, 1)
+        time_end = datetime.datetime(2014, 12, 31)
+
+        ds_tslice = dsData.sel(time=slice(time_start, time_end))
+
+        # The following line converts from DASK to numpy to supress an odd
+        # warning that doesn't influence the figure output
+        ds_tslice.SSS.values
+
+        monthly_clim_data = ds_tslice.groupby('time.month').mean('time')
+
+        # Rename the observation data for code compactness
+        dsData = monthly_clim_data.transpose('month', 'lon', 'lat')
+        obsFieldName = 'SSS'
+
+        # Set appropriate figure labels for SSS
+        preIndustrial_txt = "2011-2014"
+
+        obsTitleLabel = "Observations (Aquarius, {})".format(preIndustrial_txt)
+        fileOutLabel = 'sssAquarius'
+        unitsLabel = 'PSU'
+
     ds = xr.open_mfdataset(
         infiles,
         preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset,
diff --git a/mpas_analysis/ocean/variable_stream_map.py b/mpas_analysis/ocean/variable_stream_map.py
index 207944cb1..1bccabb7a 100644
--- a/mpas_analysis/ocean/variable_stream_map.py
+++ b/mpas_analysis/ocean/variable_stream_map.py
@@ -20,7 +20,7 @@
 oceanVariableMap['avgLayerTemperature'] = \
     ['time_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature',
      'time_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature_1',
-    'timeMonthly_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature']
+     'timeMonthly_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature']
 oceanVariableMap['sumLayerMaskValue'] = \
     ['time_avg_avgValueWithinOceanLayerRegion_sumLayerMaskValue',
      'time_avg_avgValueWithinOceanLayerRegion_sumLayerMaskValue_1',
@@ -44,3 +44,8 @@
     ['time_avg_activeTracers_temperature',
      'time_avg_activeTracers_temperature_1',
      'timeMonthly_avg_activeTracers_temperature']
+
+oceanVariableMap['sss'] = \
+    ['time_avg_activeTracers_salinity',
+     'time_avg_activeTracers_salinity_1',
+     'timeMonthly_avg_activeTracers_salinity']
diff --git a/run_analysis.py b/run_analysis.py
index 864d77327..493bdfbbc 100644
--- a/run_analysis.py
+++ b/run_analysis.py
@@ -137,6 +137,14 @@ def analysis(config):  # {{{
         ocn_modelvsobs(config, 'mld', streamMap=oceanStreamMap,
                        variableMap=oceanVariableMap)
 
+    if config.getboolean('sss_modelvsobs', 'generate'):
+        print ""
+        print "Plotting 2-d maps of SSS climatologies..."
+        from mpas_analysis.ocean.ocean_modelvsobs import ocn_modelvsobs
+        ocn_modelvsobs(config, 'sss', streamMap=oceanStreamMap,
+                       variableMap=oceanVariableMap)
+
+
     # GENERATE SEA-ICE DIAGNOSTICS
     if config.getboolean('seaice_timeseries', 'generate'):
         print ""