diff --git a/plotting/operationalPlots.py b/plotting/operationalPlots.py
index a001b53..3ea5063 100644
--- a/plotting/operationalPlots.py
+++ b/plotting/operationalPlots.py
@@ -533,36 +533,36 @@ def plotWaveProfile(x, waveHs, bathyToPlot, fname):
     plt.close()
 
 
-# these are all the ones that were formerly in CSHORE_plotLib
-def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
-    """
-    This script basically just compares two time series, under
-        the assmption that one is from the model and one a set of observations
+def adjust_plot_ticks(p_dict):
+    """ This scripts adjusts the axis scale factor of the time series plot in obs_v_mod_ts.
 
-    :param  file_path: this is the full file-path (string) to the location where the plot will be saved
-    :param p_dict: has 6 keys to it.
-        (1) a vector of datetimes ('time')
-        (2) vector of observations ('obs')
-        (3) vector of model data ('model')
-        (4) variable name (string) ('var_name')
-        (5) variable units (string!!) ('units') -> this will be put inside a tex math environment!!!!
-        (6) plot title (string) ('p_title')
-    :return: a model vs. observation time-series plot'
-        the dictionary of the statistics calculated
+           Args:
 
-    """
-    # this function plots observed data vs. model data for time-series data and computes stats for it.
+            p_dict: has 6 keys to it.
+              'time':  a vector of datetimes
 
-    assert len(p_dict['time']) == len(p_dict['obs']) == len(p_dict['model']), "Your time, model, and observation arrays are not all the same length!"
-    assert sum([isinstance(p_dict['time'][ss], DT.datetime) for ss in range(0, len(p_dict['time']))]) == len(p_dict['time']), 'Your times input must be an array of datetimes!'
+              'obs':  vector of observations
+
+              'model': vector of model data
+
+              'var_name' (str): variable name
+
+              'units' (str): variable units -> this will be put inside a tex math environment!!!!
+
+              'p_title' (str): plot title"""
+
+    assert len(p_dict['time']) == len(p_dict['obs']) == len(
+        p_dict['model']), "Your time, model, and observation arrays are not all the same length!"
+    assert sum([isinstance(p_dict['time'][ss], DT.datetime) for ss in range(0, len(p_dict['time']))]) == len(
+        p_dict['time']), 'Your times input must be an array of datetimes!'
     # calculate total duration of data to pick ticks for Xaxis on time series plot
     totalDuration = p_dict['time'][-1] - p_dict['time'][0]
     if totalDuration.days > 365:  # this is a year +  of data
         # mark 7 day increments with monthly major lables
-        majorTickLocator = mdates.MonthLocator(interval=3) # every 3 months
-        minorTickLocator = mdates.AutoDateLocator() # DayLocator(7)
+        majorTickLocator = mdates.MonthLocator(interval=3)  # every 3 months
+        minorTickLocator = mdates.AutoDateLocator()  # DayLocator(7)
         xfmt = mdates.DateFormatter('%Y-%m')
-    elif totalDuration.days > 30: # thie is months of data that is not a year
+    elif totalDuration.days > 30:  # thie is months of data that is not a year
         # mark 12 hour with daily major labels
         majorTickLocator = mdates.DayLocator(1)
         minorTickLocator = mdates.HourLocator(12)
@@ -574,30 +574,36 @@ def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
         xfmt = mdates.DateFormatter('%Y-%m-%d %H:%M')
     else:
         # mark hourly with 6 hour labels major intervals
-        tickInterval = 12  # hours?
+        tickInterval = 6  # hours?
         majorTickLocator = mdates.HourLocator(interval=tickInterval)
         minorTickLocator = mdates.HourLocator(1)
-        xfmt = mdates.DateFormatter('%m/%d\n%H:%M')
-    # DLY notes 12/17/2018 - I think this tick selection section still needs work,
-    # it works fine in some cases but terrible in others
+        xfmt = mdates.DateFormatter('%Y-%m-%d %H:%M')
 
-    ####################################################################################################################
-    # Begin Plot
-    ####################################################################################################################
-    fig = plt.figure(figsize=(10, 10))
-    if 'p_title' in p_dict.keys():
-        fig.suptitle(p_dict['p_title'], fontsize=18, fontweight='bold', verticalalignment='top')
+    return xfmt, minorTickLocator, majorTickLocator
 
-    # time series
-    ax1 = plt.subplot2grid((2, 2), (0, 0), colspan=2)
+
+def determine_axis_scale_factor(p_dict):
+    """ This scripts determines the axis scale factor of the time series plot in obs_v_mod_ts.
+
+       Args:
+
+        p_dict: has 6 keys to it.
+          'time':  a vector of datetimes
+
+          'obs':  vector of observations
+
+          'model': vector of model data
+
+          'var_name' (str): variable name
+
+          'units' (str): variable units -> this will be put inside a tex math environment!!!!
+
+          'p_title' (str): plot title"""
+
+    # determine axis scale factor
     min_val = np.nanmin([np.nanmin(p_dict['obs']), np.nanmin(p_dict['model'])])
     max_val = np.nanmax([np.nanmax(p_dict['obs']), np.nanmax(p_dict['model'])])
-    if min_val < 0 and max_val > 0:
-        ax1.plot(p_dict['time'], np.zeros(len(p_dict['time'])), 'k--')
-    ax1.plot(p_dict['time'], p_dict['obs'], 'r.', label='Observed')
-    ax1.plot(p_dict['time'], p_dict['model'], 'b.', label='Model')
-    ax1.set_ylabel('%s [$%s$]' % (p_dict['var_name'], p_dict['units']), fontsize=16)
-    # determine axis scale factor
+
     if min_val >= 0:
         sf1 = 0.9
     else:
@@ -606,24 +612,289 @@ def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
         sf2 = 1.1
     else:
         sf2 = 0.9
+
+    return min_val, max_val, sf1, sf2
+
+
+def plot_string_message(p_dict, stats_dict):
+    """ This script creates the statistics string used for the one:one plot in obs_v_mod_ts
+
+       Args:
+
+        p_dict: has 6 keys to it.
+          'time':  a vector of datetimes
+
+          'obs':  vector of observations
+
+          'model': vector of model data
+
+          'var_name' (str): variable name
+
+          'units' (str): variable units -> this will be put inside a tex math environment!!!!
+
+          'p_title' (str): plot title"""
+
+
+    header_str = '%s Comparison \nModel to Observations:' % (p_dict['var_name'])
+    m_mean_str = '\n Model Mean $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['m_mean']), p_dict['units'])
+    o_mean_str = '\n Observation Mean $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['o_mean']), p_dict['units'])
+    bias_str = '\n Bias $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['bias']), p_dict['units'])
+    RMSE_str = '\n RMSE $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['RMSE']), p_dict['units'])
+    SI_str = '\n Similarity Index $=%s$' % ("{0:.2f}".format(stats_dict['scatterIndex']))
+    sym_slp_str = '\n Symmetric Slope $=%s$' % ("{0:.2f}".format(stats_dict['symSlope']))
+    corr_coef_str = '\n Correlation Coefficient $=%s$' % ("{0:.2f}".format(stats_dict['corr']))
+    RMSE_Norm_str = '\n %%RMSE $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['RMSEnorm']), p_dict['units'])
+    num_String = '\n Number of samples $= %s$' % len(stats_dict['residuals'])
+    plot_str = m_mean_str + o_mean_str + bias_str + RMSE_str + SI_str + sym_slp_str + corr_coef_str + RMSE_Norm_str + num_String
+
+    return plot_str, header_str
+
+
+def create_obs_dictionary(d1, d2, stations):
+    """This script returns a structured dictionary of several stations of interests.
+
+    Args:
+     d1 = start time
+     d2 = end time
+     stations = list of stations of interest
+
+    Returns:
+      obs_dict = has 4 keys; each being a dictionary of the stations of interests.
+      Each station key dictionaries will have 4 keys:
+          time = Time
+          Hs = Significant wave height
+          WaveDm = Mean wave direction
+          WaveTm = Mean sea surface wave period"""
+
+    go = getDataFRF.getObs(d1, d2)  # call get observations class (from getdatatestbed's getDataFRF)
+
+    # create field dictionary
+    obs_dict = {}
+
+    for i, sta in enumerate(stations):
+
+        rawspec = go.getWaveSpec(gaugenumber=sta)  # call get wave spectra function
+
+        if rawspec is not None:  # check if station has data for the provided date
+
+            obs_dict[sta] = {'time': rawspec['time'],  # time
+                             'Hs': rawspec['Hs'],  # Significant wave height
+                             'waveTm': rawspec['Tm'],  # Mean sea surface wave period
+                             'waveDm': rawspec['waveDm']}  # Mean wave direction
+        else:
+            print("{}'s rawspec is NoneType".format(sta))
+
+    return obs_dict
+
+
+def create_model_dictionary(d1, d2, stations):
+    """This script returns a structured dictionary of models and their results at stations of interests.
+
+    Args:
+     d1 = start time
+     d2 = end time
+     stations = list of stations of interest
+
+    Returns:
+      model_dict = has 4 keys; 2 being the stations and model least, and the final two each being a dictionary
+      of the models CMS and STWAVE.
+      Each model key will have station dictionaries containing:
+          time = Time
+          Hs = Significant wave height
+          WaveDm = Mean wave direction
+          WaveTm = Mean sea surface wave period"""
+
+    models = ['STWAVE', 'CMS']
+
+    gm = getDataFRF.getDataTestBed(d1, d2)  # call get data testbed class (from getdatatestbed's getDataFRF)
+
+    # create model dictionary
+    model_dict = {}
+
+    for j, mod in enumerate(models):
+
+        model_dict[mod] = {}
+
+        for i, sta in enumerate(stations):
+
+            rawspec = gm.getWaveSpecModel(prefix='HP', model=mod, gaugenumber=sta)  # call get wave spectra function
+
+            if rawspec is not None:  # check if model's station has data for the provided date
+
+                model_dict[mod][sta] = {'time': rawspec['time'],  # time
+                                        'Hs': rawspec['Hs'],  # Significant wave height
+                                        'waveTm': rawspec['waveTm'],  # Mean sea surface wave period
+                                        'waveDm': rawspec['waveDm']}  # Mean wave direction
+
+
+            else:
+                print("{}'s rawspec is NoneType".format(sta))
+
+    return model_dict
+
+
+def statistics_dictionary(observation, model, var):
+    """ This script calls statsBryant (both non-directioanl and directional versions) function and creates stats_dict.
+
+     Args:
+
+        'observation':  vector of observationervations
+
+        'model': vector of model data
+
+        'var': name of variable
+
+    Returns:
+        stats_dict: has 14 keys to it:
+
+            'bias',
+
+            'RMSEdemeaned',
+
+            'RMSE',
+
+            'RMSEnorm',
+
+            'scatterIndex',
+
+            'symSlope',
+
+            'corr',
+
+            'r2',
+
+            'PscoreWilmont',
+
+            'PscoreIMEDS',
+
+            'residuals',
+
+            'm_mean',
+
+            'o_mean'"""
+
+    stats_dict = {}
+
+    if isinstance(observation, np.ma.masked_array) and ~observation.mask.any():
+        observation = np.array(observation)
+    if isinstance(model, np.ma.masked_array) and ~model.mask.any():
+        model = np.array(model)
+
+    if 'D' in var:  # check if var is directional, to call circular stats function
+        stats_dict = statsBryantCirc(observation, model)
+    else:  # check if var is non-directional
+        stats_dict = statsBryant(observation, model)
+
+    stats_dict['m_mean'] = np.nanmean(model)  # compute model mean
+    stats_dict['o_mean'] = np.nanmean(observation)  # compute observation mean
+
+    del stats_dict['meta']  # remove meta key from stats_dict
+
+    return stats_dict
+
+
+def statistics_dictionary_p_dict(p_dict):
+    """ This script calls statsBryant function and creates stats_dict (used in obs_v_mod_ts,
+     obs_v_mod_bathy, and obs_v_mod_bathy_TN)
+
+     Args:
+
+      p_dict: has 6 keys to it.
+        'time':  a vector of datetimes
+
+        'obs':  vector of observations
+
+        'model': vector of model data
+
+        'var_name' (str): variable name
+
+        'units' (str): variable units -> this will be put inside a tex math environment!!!!
+
+        'p_title' (str): plot title"""
+
+    stats_dict = {}
+    if isinstance(p_dict['obs'], np.ma.masked_array) and ~p_dict['obs'].mask.any():
+        p_dict['obs'] = np.array(p_dict['obs'])
+    stats_dict = statsBryant(p_dict['obs'], p_dict['model'])
+    stats_dict['m_mean'] = np.nanmean(p_dict['model'])
+    stats_dict['o_mean'] = np.nanmean(p_dict['obs'])
+    del stats_dict['meta']  # remove meta key from stats_dict
+
+    return stats_dict
+
+
+def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
+    """This script basically just compares two time series, under
+        the assmption that one is from the model and one a set of observations
+
+    Args:
+      file_path: this is the full file-path (string) to the location where the plot will be saved
+      p_dict: has 6 keys to it.
+        'time':  a vector of datetimes
+
+        'obs':  vector of observations
+
+        'model': vector of model data
+
+        'var_name' (str): variable name
+
+        'units' (str): variable units -> this will be put inside a tex math environment!!!!
+
+        'p_title' (str): plot title
+
+      ofname: output file name
+      logo_path: path to a small logo to put at the bottom of the figure (Default value = 'ArchiveFolder/CHL_logo.png')
+
+    Returns:
+      a model vs. observation time-series plot'
+      the dictionary of the statistics calculated
+
+    """
+    # this function plots observed data vs. model data for time-series data and computes stats for it.
+
+    ####################################################################################################################
+    # Begin Plot
+    ####################################################################################################################
+    fig = plt.figure(figsize=(10, 10))
+    fig.suptitle(p_dict['p_title'], fontsize=18, fontweight='bold', verticalalignment='top')
+
+    # time series
+    ax1 = plt.subplot2grid((2, 2), (0, 0), colspan=2)
+
+    min_val, max_val, sf1, sf2 = determine_axis_scale_factor(p_dict)
+    xfmt, minorTickLocator, majorTickLocator = adjust_plot_ticks(p_dict)
+
+    if min_val < 0 and max_val > 0:
+        base_date = min(p_dict['time']) - DT.timedelta(
+            seconds=0.5 * (p_dict['time'][1] - p_dict['time'][0]).total_seconds())
+        base_times = np.array(
+            [base_date + DT.timedelta(seconds=n * (p_dict['time'][1] - p_dict['time'][0]).total_seconds()) for n in
+             range(0, len(p_dict['time']) + 1)])
+        ax1.plot(base_times, np.zeros(len(base_times)), 'k--')
+
+    plt.grid()
+    ax1.scatter(p_dict['time'], p_dict['obs'], s=75, c='r', marker='o', label='Observed')
+    ax1.scatter(p_dict['time'], p_dict['model'], s=75, c='b', marker='o', label='Model')
+    ax1.set_ylabel('%s [$%s$]' % (p_dict['var_name'], p_dict['units']), fontsize=16)
+
     ax1.set_ylim([sf1 * min_val, sf2 * max_val])
-    ax1.set_xlim([min(p_dict['time']) - DT.timedelta(seconds=0.5 * (p_dict['time'][1] - p_dict['time'][0]).total_seconds()),
-                  max(p_dict['time']) + DT.timedelta(seconds=0.5 * (p_dict['time'][1] - p_dict['time'][0]).total_seconds())])
+    ax1.set_xlim(
+        [min(p_dict['time']) - DT.timedelta(seconds=0.5 * (p_dict['time'][1] - p_dict['time'][0]).total_seconds()),
+         max(p_dict['time']) + DT.timedelta(seconds=0.5 * (p_dict['time'][1] - p_dict['time'][0]).total_seconds())])
 
     # this is what you change for time-series x-axis ticks!!!!!
-    #
     # ax1.xaxis.set_major_locator(majorTickLocator)
     # ax1.xaxis.set_minor_locator(minorTickLocator)
     # ax1.xaxis.set_major_formatter(xfmt)
+
     for tick in ax1.xaxis.get_major_ticks():
         tick.label.set_fontsize(14)
     for tick in ax1.yaxis.get_major_ticks():
         tick.label.set_fontsize(14)
 
-    ax1.minorticks_off()
     ax1.tick_params(labelsize=14)
-    plt.legend(bbox_to_anchor=(0., 1.02, 1., 0.102), loc=3, ncol=3, borderaxespad=0., fontsize=14)
-    fig.autofmt_xdate()
+    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=3, borderaxespad=0., fontsize=14)
+
     # Now working on the 1-1 comparison subplot
     one_one = np.linspace(min_val - 0.05 * (max_val - min_val), max_val + 0.05 * (max_val - min_val), 100)
     ax2 = plt.subplot2grid((2, 2), (1, 0), colspan=1)
@@ -631,7 +902,7 @@ def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
     if min_val < 0 and max_val > 0:
         ax2.plot(one_one, np.zeros(len(one_one)), 'k--')
         ax2.plot(np.zeros(len(one_one)), one_one, 'k--')
-    ax2.plot(p_dict['obs'], p_dict['model'], 'r*')
+    ax2.scatter(p_dict['obs'], p_dict['model'], s=75, c='r', marker='*')
     ax2.set_xlabel('Observed %s [$%s$]' % (p_dict['var_name'], p_dict['units']), fontsize=16)
     ax2.set_ylabel('Model %s [$%s$]' % (p_dict['var_name'], p_dict['units']), fontsize=16)
     ax2.set_xlim([min_val - 0.025 * (max_val - min_val), max_val + 0.025 * (max_val - min_val)])
@@ -643,45 +914,220 @@ def obs_V_mod_TS(ofname, p_dict, logo_path='ArchiveFolder/CHL_logo.png'):
     ax2.tick_params(labelsize=14)
     plt.legend(loc=0, ncol=1, borderaxespad=0.5, fontsize=14)
 
-    # stats and stats text
-    stats_dict = statsBryant(p_dict['obs'], p_dict['model'])
-    stats_dict['m_mean'] = np.nanmean(p_dict['model'])
-    stats_dict['o_mean'] = np.nanmean(p_dict['obs'])
+    stats_dict = statistics_dictionary_p_dict(p_dict)
+    plot_str, header_str = plot_string_message(p_dict, stats_dict)
 
-    header_str = '%s Comparison \nModel to Observations:' % (p_dict['var_name'])
-    m_mean_str = '\n Model Mean $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['m_mean']), p_dict['units'])
-    o_mean_str = '\n Observation Mean $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['o_mean']), p_dict['units'])
-    bias_str = '\n Bias $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['bias']), p_dict['units'])
-    RMSE_str = '\n RMSE $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['RMSE']), p_dict['units'])
-    SI_str = '\n Similarity Index $=%s$' % ("{0:.2f}".format(stats_dict['scatterIndex']))
-    sym_slp_str = '\n Symmetric Slope $=%s$' % ("{0:.2f}".format(stats_dict['symSlope']))
-    corr_coef_str = '\n Correlation Coefficient $=%s$' % ("{0:.2f}".format(stats_dict['corr']))
-    RMSE_Norm_str = '\n %%RMSE $=%s$ $(%s)$' % ("{0:.2f}".format(stats_dict['RMSEnorm']), p_dict['units'])
-
-    num_String = '\n Number of samples $= %s$' %len(stats_dict['residuals'])
-    plot_str = m_mean_str + o_mean_str + bias_str + RMSE_str + RMSE_Norm_str + SI_str + sym_slp_str + corr_coef_str + num_String
     ax3 = plt.subplot2grid((2, 2), (1, 1), colspan=1)
+    ax3.axis('off')
     ax4 = ax3.twinx()
     ax3.axis('off')
-    ax4.axis('off')
+
     try:
-        CHL_logo = image.imread(logo_path)
+        ax4.axis('off')
+        CHL_logo = image.imread(os.path.join(logo_path))
         ax4 = fig.add_axes([0.78, 0.02, 0.20, 0.20], anchor='SE', zorder=-1)
         ax4.imshow(CHL_logo)
         ax4.axis('off')
     except:
         print('Plot generated sans CHL Logo!')
 
+    ax3.axis('off')
     ax3.text(0.01, 0.99, header_str, verticalalignment='top', horizontalalignment='left', color='black', fontsize=18,
              fontweight='bold')
     ax3.text(0.01, 0.90, plot_str, verticalalignment='top', horizontalalignment='left', color='black', fontsize=16)
 
     fig.subplots_adjust(wspace=0.4, hspace=0.1)
-    # fig.tight_layout(pad=1, h_pad=2.5, w_pad=1, rect=[0.0, 0.0, 1.0, 0.925])
-    fig.savefig(ofname, dpi=300)
+    fig.tight_layout(pad=1, h_pad=2.5, w_pad=1, rect=[0.0, 0.0, 1.0, 0.925])
+    fig.savefig(ofname, dpi=80)
     plt.close()
+
     return stats_dict
 
+
+def obs_v_mod_subplot(ofname,obs_dict, model_dict):
+    """ This script compares model's time series to observational data. Each row is a different station;
+    whereas each column is a different variable (for example: Hs, Tm, Dir, etc...).
+
+            Args:
+                ofname: output file name
+                obs_dict: dictionary with stations observations
+                model_dict: dictionary with model results at station locations"""
+
+    # figure size
+    width = 14  # figure width
+    height = 10  # figure height
+
+    dontPlot = ['time', 'units']  # variables we don't want to plot
+    colors = ['r', 'b', 'g', 'c','magenta','orange']  # line colors for models
+    obslen = len(obs_dict.keys())  # number of obervations (rows in plot)
+
+    varn = []  # variable list (Hs, waveTm, etc...)
+    if not varn:
+        varn = list(list(obs_dict.keys()))[0]
+
+    # create ylabel and unit lists
+    ylabel = []
+    for temp in obs_dict[varn].keys():
+        if temp not in dontPlot:
+            ylabel.append(temp)
+        elif temp == 'units':
+            units = list(obs_dict[varn][temp].values()) # ylabel unit list
+    varlen = len(ylabel)
+
+    ## create time series subplot figure
+    f, ax = plt.subplots(obslen, varlen, sharex=True, figsize=(width, height))
+
+    ## plot observation time series
+    for ii, obs in enumerate(obs_dict.keys()):
+        for jj, var in enumerate(obs_dict[obs].keys()):
+            if var not in dontPlot:
+                ax[ii, jj - len(dontPlot)].plot(obs_dict[obs]['time'], obs_dict[obs][var], 'k.')
+
+    ## plot Model time series
+    for jj, mod in enumerate(model_dict.keys()):
+        for ii, obs in enumerate(obs_dict.keys()):
+            if obs in model_dict[mod].keys():
+                for kk, var in enumerate(obs_dict[obs]):
+                    if var not in dontPlot:
+                        ax[ii, kk - len(dontPlot)].plot(model_dict[mod][obs]['time'], model_dict[mod][obs][var],
+                                                        colors[jj])
+
+    ## Add station label per row
+    station_list = []  # create station name list
+    for temp in obs_dict.keys():
+        station_list.append(temp)
+
+    for i in range(obslen):
+        for j in range(varlen):
+            ax[i, j].text(0.70, 0.10, station_list[i],
+                          transform=ax[i, j].transAxes,
+                          horizontalalignment='left',
+                          verticalalignment='top',
+                          weight="bold",
+                          fontsize=12)
+
+            ## Add Xlabel and grid
+    for axs in ax.flat:
+        axs.set_xlabel('Date', fontsize=16)
+        axs.grid(True)
+
+    ## Add legend
+    leg = ['Obs']
+    for temp in model_dict.keys():
+        leg.append(temp)
+
+    f.autofmt_xdate()
+    f.legend([11, 12, 13],
+             labels=leg,
+             loc='upper center',
+             ncol=3,
+             bbox_to_anchor=(0.5, 1.05),
+             fontsize=16)
+
+    ## Add Ylabels
+    for kk, label in enumerate(ylabel):
+        x_coord = 1 / varlen * kk
+        plt.figtext(x_coord, 0.5, f"{label} ({units[kk]})", fontsize=16, rotation='vertical')
+
+    ## save figure
+    f.tight_layout()
+    plt.savefig(ofname, bbox_inches='tight', dpi=80)
+    plt.close()
+
+
+def obs_v_mod_stats_subplot(ofname,stats_dict):
+    """ This script compares bar charts of stats_dict.
+
+        Args:
+            ofname: output file name
+            stats_dict: has keys for each model compared to field data stations. Each variable has 14 stats keys:
+
+                'bias',
+
+                'RMSEdemeaned',
+
+                'RMSE',
+
+                'RMSEnorm',
+
+                'scatterIndex',
+
+                'symSlope',
+
+                'corr',
+
+                'r2',
+
+                'PscoreWilmont',
+
+                'PscoreIMEDS',
+
+                'residuals',
+
+                'm_mean',
+
+                'o_mean'"""
+
+    # create model, station, y-label and x-label lists
+    y_labels = []
+    if not y_labels:
+        model = list(stats_dict.keys())  # model list
+        station_list = list(stats_dict[model[0]].keys())  # field station list
+        y_labels = list(stats_dict[model[0]][station_list[0]].keys())  # y label list (variables)
+        x_labels = list(stats_dict[model[0]][station_list[0]][y_labels[0]].keys())  # x label list (stats)
+
+    varlen = len(y_labels)  # number of variables (figure column)
+    obslen = len(station_list)  # number of field stations (figure row)
+    statlen = len(x_labels)  # number of stats variables
+
+    # figure specs
+    width = 14  # figure width
+    height = 8  # figure height
+    colors = ['r', 'b', 'g', 'c', 'm', 'y']  # bar colors for models
+    ind = np.arange(statlen)  # x axis tick array
+    wid = 1 / (len(model) + 1)  # bar width
+
+    # start figure
+    f, ax = plt.subplots(obslen, varlen, sharex=True, figsize=(width, height))
+    for ii, mod in enumerate(stats_dict.keys()):
+        for jj, obs in enumerate(stats_dict[mod].keys()):
+            for kk, var in enumerate(stats_dict[mod][obs].keys()):
+                for ll, stat in enumerate(stats_dict[mod][obs][var].keys()):
+                    ax[jj, kk].bar(ind[ll] + wid * ii, stats_dict[mod][obs][var][stat], color=colors[ii],
+                                   align='center', width=wid)
+
+    # set x labels and field station labels per row
+    for i in range(obslen):
+        for j in range(varlen):
+            ax[i, j].grid(which='major', axis='y', linestyle='--') # set grid
+            ax[i, j].set_xticks(ind + wid)  # set x tick locations
+            ax[i, j].set_xticklabels(labels=x_labels, rotation=45, horizontalalignment='right')  # set x tick labels
+            ax[i, j].text(0.05, 0.90, station_list[i],  # set field station label per row
+                          transform=ax[i, j].transAxes,
+                          horizontalalignment='left',
+                          verticalalignment='top',
+                          weight="bold",
+                          fontsize=12)
+
+    # set y labels per column
+    for kk, label in enumerate(y_labels):
+        x_coord = 1 / varlen * kk
+        plt.figtext(x_coord, 0.5, label, fontsize=16, rotation='vertical')
+
+    # set legend at top of plot
+    leg = f.legend(labels=model,
+                   loc='upper center',
+                   ncol=len(model),
+                   bbox_to_anchor=(0.5, 1.05),
+                   fontsize=16)
+    for i in range(len(model)):
+        leg.legendHandles[i].set_color(colors[i])
+
+    # save figure
+    f.tight_layout()
+    plt.savefig(ofname, bbox_inches='tight', dpi=80)
+    plt.close()
+
 def bc_plot(ofname, p_dict):
     """
     This is the script to plot some information about the boundary conditions that were put into the CSHORE infile..