diff --git a/jobs/JREGIONAL_PLOT_ALLVARS b/jobs/JREGIONAL_PLOT_ALLVARS
new file mode 100755
index 0000000000..29f4371a4e
--- /dev/null
+++ b/jobs/JREGIONAL_PLOT_ALLVARS
@@ -0,0 +1,112 @@
+#!/bin/bash
+
+#
+#-----------------------------------------------------------------------
+#
+# Source the variable definitions file and the bash utility functions.
+#
+#-----------------------------------------------------------------------
+#
+. $USHdir/source_util_funcs.sh
+source_config_for_task "task_plot_allvars|task_run_fcst" ${GLOBAL_VAR_DEFNS_FP}
+. $USHdir/job_preamble.sh
+#
+#-----------------------------------------------------------------------
+#
+# Save current shell options (in a global array).  Then set new options
+# for this script/function.
+#
+#-----------------------------------------------------------------------
+#
+{ save_shell_opts; . $USHdir/preamble.sh; } > /dev/null 2>&1
+#
+#-----------------------------------------------------------------------
+#
+# Get the full path to the file in which this script/function is located 
+# (scrfunc_fp), the name of that file (scrfunc_fn), and the directory in
+# which the file is located (scrfunc_dir).
+#
+#-----------------------------------------------------------------------
+#
+scrfunc_fp=$( $READLINK -f "${BASH_SOURCE[0]}" )
+scrfunc_fn=$( basename "${scrfunc_fp}" )
+scrfunc_dir=$( dirname "${scrfunc_fp}" )
+#
+#-----------------------------------------------------------------------
+#
+# Print message indicating entry into script.
+#
+#-----------------------------------------------------------------------
+#
+print_info_msg "
+========================================================================
+Entering script:  \"${scrfunc_fn}\"
+In directory:     \"${scrfunc_dir}\"
+
+This is the J-job script for the plotting task
+========================================================================"
+#
+#-----------------------------------------------------------------------
+#
+# Set grid name and COMOUT locations.
+#
+#-----------------------------------------------------------------------
+#
+GRID_NAME=${PREDEF_GRID_NAME:-$EXPT_SUBDIR}
+if [ "${RUN_ENVIR}" != "nco" ]; then
+  export COMOUT="${COMOUT}${SLASH_ENSMEM_SUBDIR}/postprd"
+fi
+# if COMOUT_REF is templatized with $PDY and $cyc, evaluate it
+# This is necessary to compare multiple cycles correctly
+COMOUT_REF=$(eval echo ${COMOUT_REF})
+#
+#-----------------------------------------------------------------------
+#
+# Call the ex-script for this J-job and pass to it the necessary variables.
+#
+#-----------------------------------------------------------------------
+#
+
+# plot all variables
+$SCRIPTSdir/exregional_plot_allvars.py \
+           --cycle ${CDATE} \
+           --start ${PLOT_FCST_START} \
+           --end ${PLOT_FCST_END:-$FCST_LEN_HRS} \
+           --inc ${PLOT_FCST_INC:-1} \
+           --comout ${COMOUT} \
+           --cartopy-dir ${FIXshp} \
+           --domain ${GRID_NAME} || \
+print_err_msg_exit "\
+Call to ex-script corresponding to J-job \"${scrfunc_fn}\" failed."
+
+# plot all variables diff with baseline
+if [ ! -z $COMOUT_REF ]; then
+  $SCRIPTSdir/exregional_plot_allvars_diff.py \
+           --cycle ${CDATE} \
+           --start ${PLOT_FCST_START} \
+           --end ${PLOT_FCST_END:-$FCST_LEN_HRS} \
+           --inc ${PLOT_FCST_INC:-1} \
+           --comout-1 ${COMOUT} \
+           --comout-2 ${COMOUT_REF} \
+           --cartopy-dir ${FIXshp} \
+           --domain ${GRID_NAME} || \
+  print_err_msg_exit "\
+  Call to ex-script corresponding to J-job \"${scrfunc_fn}\" failed."
+fi
+#
+#-----------------------------------------------------------------------
+#
+# Run job postamble.
+#
+#-----------------------------------------------------------------------
+#
+job_postamble
+#
+#-----------------------------------------------------------------------
+#
+# Restore the shell options saved at the beginning of this script/function.
+#
+#-----------------------------------------------------------------------
+#
+{ restore_shell_opts; } > /dev/null 2>&1
+
diff --git a/modulefiles/tasks/cheyenne/plot_allvars.local.lua b/modulefiles/tasks/cheyenne/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/cheyenne/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/gaea/plot_allvars.local.lua b/modulefiles/tasks/gaea/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/gaea/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/hera/plot_allvars.local.lua b/modulefiles/tasks/hera/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/hera/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/jet/plot_allvars.local.lua b/modulefiles/tasks/jet/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/jet/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/noaacloud/plot_allvars.local.lua b/modulefiles/tasks/noaacloud/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/noaacloud/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/orion/plot_allvars.local.lua b/modulefiles/tasks/orion/plot_allvars.local.lua
new file mode 100644
index 0000000000..d371c08301
--- /dev/null
+++ b/modulefiles/tasks/orion/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("miniconda_regional_workflow")
diff --git a/modulefiles/tasks/wcoss2/plot_allvars.local.lua b/modulefiles/tasks/wcoss2/plot_allvars.local.lua
new file mode 100644
index 0000000000..3370fa018b
--- /dev/null
+++ b/modulefiles/tasks/wcoss2/plot_allvars.local.lua
@@ -0,0 +1 @@
+load("python_regional_workflow")
diff --git a/parm/FV3LAM_wflow.xml b/parm/FV3LAM_wflow.xml
index b0b632fee4..c1cb9bd171 100644
--- a/parm/FV3LAM_wflow.xml
+++ b/parm/FV3LAM_wflow.xml
@@ -67,6 +67,7 @@ Workflow task names.
 <!ENTITY VX_ENSPOINT_TN            "{{ vx_enspoint_tn }}">
 <!ENTITY VX_ENSPOINT_MEAN_TN       "{{ vx_enspoint_mean_tn }}">
 <!ENTITY VX_ENSPOINT_PROB_TN       "{{ vx_enspoint_prob_tn }}">
+<!ENTITY PLOT_ALLVARS_TN           "{{ plot_allvars_tn }}">
 
 <!--
 Flags that specify whether to run the preprocessing and/or verification tasks.
@@ -755,6 +756,44 @@ the <task> tag to be identical to the ones above for other output times.
   {%- endif %}
 {%- endif %}
 
+{%- if run_task_plot_allvars %}
+<!--
+************************************************************************
+************************************************************************
+-->
+  <task name="&PLOT_ALLVARS_TN;" cycledefs="forecast" maxtries="{{ maxtries_plot_allvars }}">
+
+    &RSRV_DEFAULT;
+    <command>&LOAD_MODULES_RUN_TASK_FP; "&PLOT_ALLVARS_TN;" "&JOBSdir;/JREGIONAL_PLOT_ALLVARS"</command>
+    <nodes>{{ nnodes_plot_allvars }}:ppn={{ ppn_plot_allvars }}</nodes>
+    <walltime>{{ wtime_plot_allvars }}</walltime>
+    <nodesize>&NCORES_PER_NODE;</nodesize>
+  {%- if machine not in ["WCOSS2"] %}
+    <native>&SCHED_NATIVE_CMD;</native>
+  {%- endif %}
+    <jobname>&PLOT_ALLVARS_TN;</jobname>
+    <join>&LOGDIR;/&PLOT_ALLVARS_TN;<cyclestr>_@Y@m@d@H</cyclestr>&LOGEXT;</join>
+
+    <envar><name>GLOBAL_VAR_DEFNS_FP</name><value>&GLOBAL_VAR_DEFNS_FP;</value></envar>
+    <envar><name>USHdir</name><value>&USHdir;</value></envar>
+    <envar><name>PDY</name><value><cyclestr>@Y@m@d</cyclestr></value></envar>
+    <envar><name>cyc</name><value><cyclestr>@H</cyclestr></value></envar>
+    <envar><name>subcyc</name><value><cyclestr>@M</cyclestr></value></envar>
+    <envar><name>LOGDIR</name><value>&LOGDIR;</value></envar>
+    <envar><name>SLASH_ENSMEM_SUBDIR</name><value><cyclestr>{{ slash_ensmem_subdir }}</cyclestr></value></envar>
+    <envar><name>ENSMEM_INDX</name><value><cyclestr>#{{ ensmem_indx_name }}#</cyclestr></value></envar>
+
+    <dependency>
+        {%- if write_dopost %}
+      <taskdep task="&RUN_FCST_TN;{{ uscore_ensmem_name }}"/>
+        {%- else %}
+      <metataskdep metatask="&RUN_POST_TN;{{ uscore_ensmem_name }}"/>
+        {%- endif %}
+    </dependency>
+
+  </task>
+{%- endif %}
+
 {%- if run_task_get_obs_ccpa %}
 <!--
 ************************************************************************
diff --git a/scripts/exregional_plot_allvars.py b/scripts/exregional_plot_allvars.py
new file mode 100755
index 0000000000..a81ed9252d
--- /dev/null
+++ b/scripts/exregional_plot_allvars.py
@@ -0,0 +1,1411 @@
+#!/usr/bin/env python3
+
+################################################################################
+####  Python Script Documentation Block
+#
+# Script name:       	plot_allvars.py
+# Script description:  	Generates plots from FV3-LAM post processed grib2 output
+# 			over the CONUS
+#
+# Authors:  Ben Blake		Org: NOAA/NWS/NCEP/EMC		Date: 2020-05-07
+#           David Wright 	Org: University of Michigan
+#
+# Instructions:		Make sure all the necessary modules can be imported.
+#                       The following command line arguments are needed:
+#                       1. Cycle date/time in YYYYMMDDHH format
+#                       2. Starting forecast hour
+#                       3. Ending forecast hour
+#                       4. Forecast hour increment
+#                       5. COMOUT: Postprocessed data location
+#                       6. CARTOPY_DIR:  Base directory of cartopy shapefiles
+#                          -Shapefiles cannot be directly downloaded to NOAA
+#                            machines from the internet, so shapefiles need to
+#                            be downloaded if geopolitical boundaries are
+#                            desired on the maps.
+#                          -File structure should be:
+#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
+#                          -More information regarding files needed to setup
+#                            display maps in Cartopy, see SRW App Users' Guide
+#                       7. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
+#                          used in forecast and in constructing the names
+#                          of the post output files.
+#
+#           		To create plots for a forecast on the RRFS_CONUS_25km
+#                       grid for hours 20-24 from the 5/7 00Z cycle with
+#                       hourly output:
+#                          python plot_allvars.py 2020050700 20 24 1 \
+#                          /path/to/expt_dirs/experiment/name \
+#                          /path/to/base/cartopy/maps RRFS_CONUS_25km
+#
+#                       The variable domains in this script can be set to either
+#                         'conus' for a CONUS map or 'regional' where the map
+#                         is defined from variables in the grib2 files
+#
+################################################################################
+
+# -------------Import modules --------------------------#
+import pygrib
+import cartopy.crs as ccrs
+from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
+import cartopy.feature as cfeature
+import matplotlib
+
+matplotlib.use("Agg")
+import io
+import matplotlib.pyplot as plt
+import dateutil.relativedelta, dateutil.parser
+from PIL import Image
+from matplotlib.gridspec import GridSpec
+import numpy as np
+import time, os, sys, multiprocessing
+import multiprocessing.pool
+from scipy import ndimage
+import pyproj
+import argparse
+import cartopy
+import logging
+import warnings
+
+# --------------Define some functions ------------------#
+
+
+def ndate(cdate, hours):
+    if not isinstance(cdate, str):
+        if isinstance(cdate, int):
+            cdate = str(cdate)
+        else:
+            logging.error("ERROR: Input cdate must be string or integer!")
+            sys.exit(1)
+    if not isinstance(hours, int):
+        if isinstance(hours, str):
+            hours = int(hours)
+        else:
+            logging.error("ERROR: Input delta must be string or integer!")
+            sys.exit(1)
+
+    indate = cdate.strip()
+    hh = indate[8:10]
+    yyyy = indate[0:4]
+    mm = indate[4:6]
+    dd = indate[6:8]
+    # set date/time field
+    parseme = yyyy + " " + mm + " " + dd + " " + hh
+    datetime_cdate = dateutil.parser.parse(parseme)
+    valid = datetime_cdate + dateutil.relativedelta.relativedelta(hours=+hours)
+    vyyyy = str(valid.year)
+    vm = str(valid.month).zfill(2)
+    vd = str(valid.day).zfill(2)
+    vh = str(valid.hour).zfill(2)
+    return vyyyy + vm + vd + vh
+
+
+def clear_plotables(ax, keep_ax_lst, fig):
+    #### - step to clear off old plottables but leave the map info - ####
+    if len(keep_ax_lst) == 0:
+        logging.info(
+            "clear_plotables WARNING keep_ax_lst has length 0. Clearing ALL plottables including map info!"
+        )
+    cur_ax_children = ax.get_children()[:]
+    if len(cur_ax_children) > 0:
+        for a in cur_ax_children:
+            if a not in keep_ax_lst:
+                # if the artist isn't part of the initial set up, remove it
+                a.remove()
+
+
+def compress_and_save(filename):
+    #### - compress and save the image - ####
+    ram = io.BytesIO()
+    plt.savefig(ram, format="png", bbox_inches="tight", dpi=150)
+    ram.seek(0)
+    im = Image.open(ram)
+    im2 = im.convert("RGB")  # .convert('P', palette=Image.ADAPTIVE)
+    im2.save(filename, format="PNG")
+
+
+def cmap_t2m():
+    # Create colormap for 2-m temperature
+    # Modified version of the ncl_t2m colormap
+
+    #fmt off
+    r = np.array(
+        [255,128,0,  70, 51, 0,  255,0, 0,  51, 255,255,255,255,255,171,128,128,36,162,255]
+    )
+    g = np.array(
+        [0,  0,  0,  70, 102,162,255,92,128,185,255,214,153,102,0,  0,  0,  68, 36,162,255]
+    )
+    b = np.array(
+        [255,128,128,255,255,255,255,0, 0,  102,0,  112,0,  0,  0,  56, 0,  68, 36,162,255]
+    )
+    #fmt on
+    xsize = np.arange(np.size(r))
+    r = r / 255.0
+    g = g / 255.0
+    b = b / 255.0
+    red = []
+    green = []
+    blue = []
+    for i in range(len(xsize)):
+        xNorm = float(i) / (float(np.size(r)) - 1.0)
+        red.append([xNorm, r[i], r[i]])
+        green.append([xNorm, g[i], g[i]])
+        blue.append([xNorm, b[i], b[i]])
+    colorDict = {"red": red, "green": green, "blue": blue}
+    cmap_t2m_coltbl = matplotlib.colors.LinearSegmentedColormap(
+        "CMAP_T2M_COLTBL", colorDict
+    )
+    return cmap_t2m_coltbl
+
+
+def cmap_q2m():
+    # Create colormap for 2-m dew point temperature
+    r = np.array([255, 179, 96, 128, 0, 0, 51, 0, 0, 0, 133, 51, 70, 0, 128, 128, 180])
+    g = np.array(
+        [255, 179, 96, 128, 92, 128, 153, 155, 155, 255, 162, 102, 70, 0, 0, 0, 0]
+    )
+    b = np.array(
+        [255, 179, 96, 0, 0, 0, 102, 155, 255, 255, 255, 255, 255, 128, 255, 128, 128]
+    )
+    xsize = np.arange(np.size(r))
+    r = r / 255.0
+    g = g / 255.0
+    b = b / 255.0
+    red = []
+    green = []
+    blue = []
+    for i in range(len(xsize)):
+        xNorm = float(i) / (float(np.size(r)) - 1.0)
+        red.append([xNorm, r[i], r[i]])
+        green.append([xNorm, g[i], g[i]])
+        blue.append([xNorm, b[i], b[i]])
+    colorDict = {"red": red, "green": green, "blue": blue}
+    cmap_q2m_coltbl = matplotlib.colors.LinearSegmentedColormap(
+        "CMAP_Q2M_COLTBL", colorDict
+    )
+    cmap_q2m_coltbl.set_over(color="deeppink")
+    return cmap_q2m_coltbl
+
+
+def rotate_wind(true_lat, lov_lon, earth_lons, uin, vin, proj, inverse=False):
+    #  Rotate winds from LCC relative to earth relative (or vice-versa if inverse==true)
+    #   This routine is vectorized and *should* work on any size 2D vg and ug arrays.
+    #   Program will quit if dimensions are too large.
+    #
+    # Input args:
+    #  true_lat = True latitidue for LCC projection (single value in degrees)
+    #  lov_lon  = The LOV value from grib (e.g. - -95.0) (single value in degrees)
+    #              Grib doc says: "Lov = orientation of the grid; i.e. the east longitude value of
+    #                              the meridian which is parallel to the Y-axis (or columns of the grid)
+    #                              along which latitude increases as the Y-coordinate increases (the
+    #                              orientation longitude may or may not appear on a particular grid).
+    #
+    #  earth_lons = Earth relative longitudes (can be an array, in degrees)
+    #  uin, vin     = Input winds to rotate
+    #
+    # Returns:
+    #  uout, vout = Output, rotated winds
+    # -----------------------------------------------------------------------------------------------------
+
+    # Get size and length of input u winds, if not 2d, raise an error
+    q = np.shape(uin)
+    ndims = len(q)
+    if ndims > 2:
+        # Raise error and quit!
+        raise SystemExit("Input winds for rotation have greater than 2 dimensions!")
+    if lov_lon > 0.0:
+        lov_lon = lov_lon - 360.0
+    dtr = np.pi / 180.0  # Degrees to radians
+
+    if not isinstance(inverse, bool):
+        raise TypeError("**kwarg inverse must be of type bool.")
+
+    # Compute rotation constant which is also
+    # known as the Lambert cone constant.  In the case
+    # of a polar stereographic projection, this is one.
+    # See the following pdf for excellent documentation
+    # http://www.dtcenter.org/met/users/docs/write_ups/velocity.pdf
+    if proj.lower() == "lcc":
+        rotcon_p = np.sin(true_lat * dtr)
+    elif proj.lower() in ["stere", "spstere", "npstere"]:
+        rotcon_p = 1.0
+    else:
+        raise SystemExit(
+            "Unsupported map projection: " + proj.lower() + " for wind rotation."
+        )
+
+    angles = rotcon_p * (earth_lons - lov_lon) * dtr
+    sinx2 = np.sin(angles)
+    cosx2 = np.cos(angles)
+
+    # Steps below are elementwise products, not matrix mutliplies
+    if inverse == False:
+        # Return the earth relative winds
+        uout = cosx2 * uin + sinx2 * vin
+        vout = -sinx2 * uin + cosx2 * vin
+    elif inverse == True:
+        # Return the grid relative winds
+        uout = cosx2 * uin - sinx2 * vin
+        vout = sinx2 * uin + cosx2 * vin
+
+    return uout, vout
+
+
+def setup_logging(debug=False):
+
+    """Calls initialization functions for logging package, and sets the
+    user-defined level for logging in the script."""
+
+    level = logging.INFO
+    if debug:
+        level = logging.DEBUG
+
+    logging.basicConfig(format="%(levelname)s: %(message)s ", level=level)
+    if debug:
+        logging.info("Logging level set to DEBUG")
+
+# -------------Start of script -------------------------#
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--cycle",
+        "-c",
+        help="Cycle date/time in YYYYMMDDHH format.",
+        required=True,
+    )
+    parser.add_argument(
+        "--start",
+        "-s",
+        help="Starting forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--end",
+        "-e",
+        help="Ending forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--inc",
+        "-i",
+        help="Increment forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--comout",
+        help="Path to directory containing post-processed files.",
+        required=True,
+    )
+    parser.add_argument(
+        "--cartopy-dir",
+        help="Path to base directory of cartopy shapefiles.",
+        required=True,
+    )
+    parser.add_argument(
+        "--domain",
+        "-d",
+        help="Name of native domain used in forecast (and in constructing post file names).",
+        required=True,
+    )
+    parser.add_argument(
+        "--debug",
+        action="store_true",
+        help="Print debug messages",
+    )
+    args = parser.parse_args()
+    
+    setup_logging(args.debug)
+    
+    # Throw away python warnings (mostly depreciation.)
+    warnings.simplefilter("ignore")
+
+    # Read date/time, forecast hour, and directory paths from command line
+    ymdh = str(args.cycle)
+    ymd = ymdh[0:8]
+    year = int(ymdh[0:4])
+    month = int(ymdh[4:6])
+    day = int(ymdh[6:8])
+    hour = int(ymdh[8:10])
+    cyc = str(hour).zfill(2)
+    logging.info(f"{year}, {month}, {day}, {hour}")
+    
+    # Define the range of forecast hours to create plots for
+    start_fhr = int(args.start)
+    end_fhr = int(args.end)
+    increment_fhr = int(args.inc)
+    if (start_fhr == end_fhr) or (increment_fhr == 0):
+        fhours = [start_fhr]
+    else:
+        num = int(((end_fhr - start_fhr) / increment_fhr) + 1)
+        fhours = np.linspace(start_fhr, end_fhr, num, dtype="int")
+    logging.info(fhours)
+    
+    COMOUT = str(args.comout)
+    CARTOPY_DIR = str(args.cartopy_dir)
+    POST_OUTPUT_DOMAIN_NAME = str(args.domain).lower()
+    
+    # Loop over forecast hours
+    for fhr in fhours:
+        fhour = str(fhr).zfill(3)
+        logging.info("Working on forecast hour " + fhour)
+        itime = ymdh
+        vtime = ndate(itime, int(fhr))
+    
+        # Define the location of the input file
+        data1 = pygrib.open(
+            COMOUT
+            + "/rrfs.t"
+            + cyc
+            + "z.prslev.f"
+            + fhour
+            + "."
+            + POST_OUTPUT_DOMAIN_NAME
+            + ".grib2"
+        )
+    
+        # Get the lats and lons
+        grids = [data1]
+        lats = []
+        lons = []
+        lats_shift = []
+        lons_shift = []
+    
+        for data in grids:
+            # Unshifted grid for contours and wind barbs
+            lat, lon = data[1].latlons()
+            lats.append(lat)
+            lons.append(lon)
+    
+            # Shift grid for pcolormesh
+            lat1 = data[1]["latitudeOfFirstGridPointInDegrees"]
+            lon1 = data[1]["longitudeOfFirstGridPointInDegrees"]
+            try:
+                nx = data[1]["Nx"]
+                ny = data[1]["Ny"]
+            except:
+                nx = data[1]["Ni"]
+                ny = data[1]["Nj"]
+            dx = data[1]["DxInMetres"]
+            dy = data[1]["DyInMetres"]
+            pj = pyproj.Proj(data[1].projparams)
+            llcrnrx, llcrnry = pj(lon1, lat1)
+            llcrnrx = llcrnrx - (dx / 2.0)
+            llcrnry = llcrnry - (dy / 2.0)
+            x = llcrnrx + dx * np.arange(nx)
+            y = llcrnry + dy * np.arange(ny)
+            x, y = np.meshgrid(x, y)
+            lon, lat = pj(x, y, inverse=True)
+            lats_shift.append(lat)
+            lons_shift.append(lon)
+    
+        # Unshifted lat/lon arrays grabbed directly using latlons() method
+        lat = lats[0]
+        lon = lons[0]
+    
+        # Shifted lat/lon arrays for pcolormesh
+        lat_shift = lats_shift[0]
+        lon_shift = lons_shift[0]
+    
+        Lat0 = data1[1]["LaDInDegrees"]
+        Lon0 = data1[1]["LoVInDegrees"]
+        logging.info(Lat0)
+        logging.info(Lon0)
+    
+        # Specify plotting domains
+        # User can add domains here, just need to specify lat/lon information below
+        # (if dom == 'conus' block)
+        domains = ["conus"]  # Other option is 'regional'
+    
+        ###################################################
+        # Read in all variables and calculate differences #
+        ###################################################
+        t1a = time.perf_counter()
+    
+        # Sea level pressure
+        slp = data1.select(name="Pressure reduced to MSL")[0].values * 0.01
+        slpsmooth = ndimage.gaussian_filter(slp, 13.78)
+    
+        # 2-m temperature
+        tmp2m = data1.select(name="2 metre temperature")[0].values
+        tmp2m = (tmp2m - 273.15) * 1.8 + 32.0
+    
+        # 2-m dew point temperature
+        dew2m = data1.select(name="2 metre dewpoint temperature")[0].values
+        dew2m = (dew2m - 273.15) * 1.8 + 32.0
+    
+        # 10-m wind speed
+        uwind = data1.select(name="10 metre U wind component")[0].values * 1.94384
+        vwind = data1.select(name="10 metre V wind component")[0].values * 1.94384
+        # Rotate winds from grid relative to Earth relative
+        uwind, vwind = rotate_wind(Lat0, Lon0, lon, uwind, vwind, "lcc", inverse=False)
+        wspd10m = np.sqrt(uwind**2 + vwind**2)
+    
+        # Surface-based CAPE
+        cape = data1.select(
+            name="Convective available potential energy", typeOfLevel="surface"
+        )[0].values
+    
+        # Surface-based CIN
+        cin = data1.select(name="Convective inhibition", typeOfLevel="surface")[0].values
+    
+        # 500 mb height, wind, vorticity
+        try:
+            z500 = data1.select(name="Geopotential Height", level=500)[0].values * 0.1
+            z500 = ndimage.gaussian_filter(z500, 6.89)
+            vort500 = data1.select(name="Absolute vorticity", level=500)[0].values * 100000
+            vort500 = ndimage.gaussian_filter(vort500, 1.7225)
+            vort500[vort500 > 1000] = 0  # Mask out undefined values on domain edge
+            u500 = data1.select(name="U component of wind", level=500)[0].values * 1.94384
+            v500 = data1.select(name="V component of wind", level=500)[0].values * 1.94384
+            # Rotate winds from grid relative to Earth relative
+            u500, v500 = rotate_wind(Lat0, Lon0, lon, u500, v500, "lcc", inverse=False)
+        except:
+            u500 = None
+            v500 = None
+    
+        # 250 mb winds
+        u250 = data1.select(name="U component of wind", level=250)[0].values * 1.94384
+        v250 = data1.select(name="V component of wind", level=250)[0].values * 1.94384
+        # Rotate winds from grid relative to Earth relative
+        u250, v250 = rotate_wind(Lat0, Lon0, lon, u250, v250, "lcc", inverse=False)
+        wspd250 = np.sqrt(u250**2 + v250**2)
+    
+        # Total precipitation
+        qpf = (
+            data1.select(name="Total Precipitation", lengthOfTimeRange=fhr)[0].values
+            * 0.0393701
+        )
+    
+        # Composite reflectivity
+        refc = data1.select(name="Maximum/Composite radar reflectivity")[0].values
+    
+        if fhr > 0:
+            # Max/Min Hourly 2-5 km Updraft Helicity
+            maxuh25 = data1.select(
+                stepType="max", parameterName="199", topLevel=5000, bottomLevel=2000
+            )[0].values
+            minuh25 = data1.select(
+                stepType="min", parameterName="200", topLevel=5000, bottomLevel=2000
+            )[0].values
+            maxuh25[maxuh25 < 10] = 0
+            minuh25[minuh25 > -10] = 0
+            uh25 = maxuh25 + minuh25
+    
+        t2a = time.perf_counter()
+        t3a = round(t2a - t1a, 3)
+        logging.info(("%.3f seconds to read all messages") % t3a)
+    
+        ########################################
+        #    START PLOTTING FOR EACH DOMAIN    #
+        ########################################
+    
+        def main():
+    
+            # Number of processes must coincide with the number of domains to plot
+            # pool = multiprocessing.Pool(len(domains))
+            # pool.map(plot_all,domains)
+    
+            # To avoid import multiprocessing recursively on MacOS etc.
+            # Anyway since we only have one domain for SRW application
+            for dom in domains:
+                plot_all(dom)
+    
+        def plot_all(dom):
+    
+            t1dom = time.perf_counter()
+    
+            # Map corners for each domain
+            if dom == "conus":
+                llcrnrlon = -120.5
+                llcrnrlat = 21.0
+                urcrnrlon = -64.5
+                urcrnrlat = 49.0
+                lat_0 = 35.4
+                lon_0 = -97.6
+                extent = [llcrnrlon - 3, urcrnrlon - 6, llcrnrlat - 1, urcrnrlat + 2]
+            elif dom == "regional":
+                llcrnrlon = np.min(lon)
+                llcrnrlat = np.min(lat)
+                urcrnrlon = np.max(lon)
+                urcrnrlat = np.max(lat)
+                lat_0 = Lat0
+                lon_0 = Lon0
+                extent = [llcrnrlon, urcrnrlon, llcrnrlat - 1, urcrnrlat]
+    
+            # create figure and axes instances
+            fig = plt.figure(figsize=(10, 10))
+            ax1 = fig.add_axes([0.1, 0.1, 0.8, 0.8])
+    
+            # Define where Cartopy Maps are located
+            cartopy.config["data_dir"] = CARTOPY_DIR
+    
+            back_res = "50m"
+            back_img = "on"
+    
+            # set up the map background with cartopy
+            myproj = ccrs.LambertConformal(
+                central_longitude=lon_0,
+                central_latitude=lat_0,
+                false_easting=0.0,
+                false_northing=0.0,
+                secant_latitudes=None,
+                standard_parallels=None,
+                globe=None,
+            )
+            ax = plt.axes(projection=myproj)
+            ax.set_extent(extent)
+    
+            fline_wd = 0.5  # line width
+            falpha = 0.3  # transparency
+    
+            # natural_earth
+            #  land=cfeature.NaturalEarthFeature('physical','land',back_res,
+            #                    edgecolor='face',facecolor=cfeature.COLORS['land'],
+            #                    alpha=falpha)
+            lakes = cfeature.NaturalEarthFeature(
+                "physical",
+                "lakes",
+                back_res,
+                edgecolor="blue",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+            coastline = cfeature.NaturalEarthFeature(
+                "physical",
+                "coastline",
+                back_res,
+                edgecolor="blue",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+            states = cfeature.NaturalEarthFeature(
+                "cultural",
+                "admin_1_states_provinces",
+                back_res,
+                edgecolor="black",
+                facecolor="none",
+                linewidth=fline_wd,
+                linestyle=":",
+                alpha=falpha,
+            )
+            borders = cfeature.NaturalEarthFeature(
+                "cultural",
+                "admin_0_countries",
+                back_res,
+                edgecolor="red",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+    
+            # All lat lons are earth relative, so setup the associated projection correct for that data
+            transform = ccrs.PlateCarree()
+    
+            # high-resolution background images
+            if back_img == "on":
+                img = plt.imread(CARTOPY_DIR + "/raster_files/NE1_50M_SR_W.tif")
+                ax.imshow(img, origin="upper", transform=transform)
+    
+            #  ax.add_feature(land)
+            ax.add_feature(lakes)
+            ax.add_feature(states)
+            ax.add_feature(borders)
+            ax.add_feature(coastline)
+    
+            # Map/figure has been set up here, save axes instances for use again later
+            keep_ax_lst = ax.get_children()[:]
+    
+            ################################
+            # Plot SLP
+            ################################
+            t1 = time.perf_counter()
+            logging.info(("Working on slp for " + dom))
+    
+            units = "mb"
+            clevs = [
+                976,
+                980,
+                984,
+                988,
+                992,
+                996,
+                1000,
+                1004,
+                1008,
+                1012,
+                1016,
+                1020,
+                1024,
+                1028,
+                1032,
+                1036,
+                1040,
+                1044,
+                1048,
+                1052,
+            ]
+            clevsdif = [-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12]
+            cm = plt.cm.Spectral_r
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs1_a = plt.pcolormesh(
+                lon_shift, lat_shift, slp, transform=transform, cmap=cm, norm=norm
+            )
+            cbar1 = plt.colorbar(
+                cs1_a, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            cs1_b = plt.contour(
+                lon_shift,
+                lat_shift,
+                slpsmooth,
+                np.arange(940, 1060, 4),
+                colors="black",
+                linewidths=1.25,
+                transform=transform,
+            )
+            plt.clabel(cs1_b, np.arange(940, 1060, 4), inline=1, fmt="%d", fontsize=8)
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM SLP ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/slp_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot slp for: " + dom) % t3)
+    
+            #################################
+            # Plot 2-m T
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on t2m for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "\xb0" "F"
+            clevs = np.linspace(-16, 134, 51)
+            cm = plt.cm.Spectral_r  # cmap_t2m()
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift, lat_shift, tmp2m, transform=transform, cmap=cm, norm=norm
+            )
+            cs_1.cmap.set_under("white")
+            cs_1.cmap.set_over("white")
+            cbar1 = plt.colorbar(
+                cs_1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=[-16, -4, 8, 20, 32, 44, 56, 68, 80, 92, 104, 116, 128],
+                extend="both",
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM 2-m Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/2mt_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 2mt for: " + dom) % t3)
+    
+            #################################
+            # Plot 2-m Dew Point
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 2mdew for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "\xb0" "F"
+            clevs = np.linspace(-5, 80, 35)
+            cm = cmap_q2m()
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift, lat_shift, dew2m, transform=transform, cmap=cm, norm=norm
+            )
+            cbar1 = plt.colorbar(
+                cs_1, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM 2-m Dew Point Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/2mdew_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 2mdew for: " + dom) % t3)
+    
+            #################################
+            # Plot 10-m WSPD
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 10mwspd for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "kts"
+            # Places a wind barb every ~180 km, optimized for CONUS domain
+            skip = round(177.28 * (dx / 1000.0) ** -0.97)
+            logging.info("skipping every " + str(skip) + " grid points to plot")
+            barblength = 4
+    
+            clevs = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
+            colorlist = [
+                "turquoise",
+                "dodgerblue",
+                "blue",
+                "#FFF68F",
+                "#E3CF57",
+                "peru",
+                "brown",
+                "crimson",
+                "red",
+                "fuchsia",
+                "DarkViolet",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift,
+                lat_shift,
+                wspd10m,
+                transform=transform,
+                cmap=cm,
+                vmin=5,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            plt.barbs(
+                lon_shift[::skip, ::skip],
+                lat_shift[::skip, ::skip],
+                uwind[::skip, ::skip],
+                vwind[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM 10-m Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/10mwind_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 10mwspd for: " + dom) % t3)
+    
+            #################################
+            # Plot Surface-Based CAPE/CIN
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on surface-based CAPE/CIN for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "J/kg"
+            clevs = [100, 250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000]
+            clevs2 = [-2000, -500, -250, -100, -25]
+            colorlist = [
+                "blue",
+                "dodgerblue",
+                "cyan",
+                "mediumspringgreen",
+                "#FAFAD2",
+                "#EEEE00",
+                "#EEC900",
+                "darkorange",
+                "crimson",
+                "darkred",
+                "darkviolet",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift,
+                lat_shift,
+                cape,
+                transform=transform,
+                cmap=cm,
+                vmin=100,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            cs_1b = plt.contourf(
+                lon_shift,
+                lat_shift,
+                cin,
+                clevs2,
+                colors="none",
+                hatches=["**", "++", "////", ".."],
+                transform=transform,
+            )
+            ax.text(
+                0.5,
+                1.05,
+                "FV3-LAM Surface-Based CAPE (shaded) and CIN (hatched) ("
+                + units
+                + ") \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/sfcape_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot surface-based CAPE/CIN for: " + dom) % t3)
+    
+            #################################
+            # Plot 500 mb HGT/WIND/VORT
+            #################################
+            if u500 is not None:
+                t1 = time.perf_counter()
+                logging.info(("Working on 500 mb Hgt/Wind/Vort for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                clear_plotables(ax, keep_ax_lst, fig)
+    
+                units = "x10${^5}$ s${^{-1}}$"
+                skip = round(177.28 * (dx / 1000.0) ** -0.97)
+                barblength = 4
+    
+                vortlevs = [16, 20, 24, 28, 32, 36, 40]
+                colorlist = ["yellow", "gold", "goldenrod", "orange", "orangered", "red"]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                norm = matplotlib.colors.BoundaryNorm(vortlevs, cm.N)
+    
+                cs1_a = plt.pcolormesh(
+                    lon_shift, lat_shift, vort500, transform=transform, cmap=cm, norm=norm
+                )
+                cs1_a.cmap.set_under("white")
+                cs1_a.cmap.set_over("darkred")
+                cbar1 = plt.colorbar(
+                    cs1_a,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    ticks=vortlevs,
+                    extend="both",
+                )
+                cbar1.set_label(units, fontsize=8)
+                cbar1.ax.tick_params(labelsize=8)
+                plt.barbs(
+                    lon_shift[::skip, ::skip],
+                    lat_shift[::skip, ::skip],
+                    u500[::skip, ::skip],
+                    v500[::skip, ::skip],
+                    length=barblength,
+                    linewidth=0.5,
+                    color="steelblue",
+                    transform=transform,
+                )
+                cs1_b = plt.contour(
+                    lon_shift,
+                    lat_shift,
+                    z500,
+                    np.arange(486, 600, 6),
+                    colors="black",
+                    linewidths=1,
+                    transform=transform,
+                )
+                plt.clabel(
+                    cs1_b, np.arange(486, 600, 6), inline_spacing=1, fmt="%d", fontsize=8
+                )
+                ax.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM 500 mb Heights (dam), Winds (kts), and $\zeta$ ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=8,
+                    transform=ax.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT + "/500_" + dom + "_f" + fhour + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot 500 mb Hgt/Wind/Vort for: " + dom) % t3)
+    
+            #################################
+            # Plot 250 mb WIND
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 250 mb WIND for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "kts"
+            skip = round(177.28 * (dx / 1000.0) ** -0.97)
+    
+            barblength = 4
+    
+            clevs = [50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150]
+            colorlist = [
+                "turquoise",
+                "deepskyblue",
+                "dodgerblue",
+                "#1874CD",
+                "blue",
+                "beige",
+                "khaki",
+                "peru",
+                "brown",
+                "crimson",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift,
+                lat_shift,
+                wspd250,
+                transform=transform,
+                cmap=cm,
+                vmin=50,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("red")
+            cbar1 = plt.colorbar(
+                cs_1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            plt.barbs(
+                lon_shift[::skip, ::skip],
+                lat_shift[::skip, ::skip],
+                u250[::skip, ::skip],
+                v250[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM 250 mb Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/250wind_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 250 mb WIND for: " + dom) % t3)
+    
+            #################################
+            # Plot Total QPF
+            #################################
+            if fhr > 0:  # Do not make total QPF plot for forecast hour 0
+                t1 = time.perf_counter()
+                logging.info(("Working on total qpf for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                clear_plotables(ax, keep_ax_lst, fig)
+    
+                units = "in"
+                clevs = [
+                    0.01,
+                    0.1,
+                    0.25,
+                    0.5,
+                    0.75,
+                    1,
+                    1.25,
+                    1.5,
+                    1.75,
+                    2,
+                    2.5,
+                    3,
+                    4,
+                    5,
+                    7,
+                    10,
+                    15,
+                    20,
+                ]
+                clevsdif = [-3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3]
+                colorlist = [
+                    "chartreuse",
+                    "limegreen",
+                    "green",
+                    "blue",
+                    "dodgerblue",
+                    "deepskyblue",
+                    "cyan",
+                    "mediumpurple",
+                    "mediumorchid",
+                    "darkmagenta",
+                    "darkred",
+                    "crimson",
+                    "orangered",
+                    "darkorange",
+                    "goldenrod",
+                    "gold",
+                    "yellow",
+                ]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+                cs_1 = plt.pcolormesh(
+                    lon_shift,
+                    lat_shift,
+                    qpf,
+                    transform=transform,
+                    cmap=cm,
+                    vmin=0.01,
+                    norm=norm,
+                )
+                cs_1.cmap.set_under("white", alpha=0.0)
+                cs_1.cmap.set_over("pink")
+                cbar1 = plt.colorbar(
+                    cs_1,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    ticks=clevs,
+                    extend="max",
+                )
+                cbar1.set_label(units, fontsize=8)
+                cbar1.ax.set_xticklabels(clevs)
+                cbar1.ax.tick_params(labelsize=8)
+                ax.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM "
+                    + fhour
+                    + "-hr Accumulated Precipitation ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=8,
+                    transform=ax.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT + "/qpf_" + dom + "_f" + fhour + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot total qpf for: " + dom) % t3)
+    
+            #################################
+            # Plot composite reflectivity
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on composite reflectivity for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            clear_plotables(ax, keep_ax_lst, fig)
+    
+            units = "dBZ"
+            clevs = np.linspace(5, 70, 14)
+            clevsdif = [20, 1000]
+            colorlist = [
+                "turquoise",
+                "dodgerblue",
+                "mediumblue",
+                "lime",
+                "limegreen",
+                "green",
+                "#EEEE00",
+                "#EEC900",
+                "darkorange",
+                "red",
+                "firebrick",
+                "darkred",
+                "fuchsia",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = plt.pcolormesh(
+                lon_shift, lat_shift, refc, transform=transform, cmap=cm, vmin=5, norm=norm
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=8)
+            cbar1.ax.tick_params(labelsize=8)
+            ax.text(
+                0.5,
+                1.03,
+                "FV3-LAM Composite Reflectivity ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT + "/refc_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot composite reflectivity for: " + dom) % t3)
+    
+            #################################
+            # Plot Max/Min Hourly 2-5 km UH
+            #################################
+            if fhr > 0:  # Do not make max/min hourly 2-5 km UH plot for forecast hour 0
+                t1 = time.perf_counter()
+                logging.info(("Working on Max/Min Hourly 2-5 km UH for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                clear_plotables(ax, keep_ax_lst, fig)
+    
+                units = "m${^2}$ s$^{-2}$"
+                clevs = [
+                    -150,
+                    -100,
+                    -75,
+                    -50,
+                    -25,
+                    -10,
+                    0,
+                    10,
+                    25,
+                    50,
+                    75,
+                    100,
+                    150,
+                    200,
+                    250,
+                    300,
+                ]
+                #   alternative colormap for just max UH if you don't want to plot the min UH too
+                #     colorlist = ['white','skyblue','mediumblue','green','orchid','firebrick','#EEC900','DarkViolet']
+                colorlist = [
+                    "blue",
+                    "#1874CD",
+                    "dodgerblue",
+                    "deepskyblue",
+                    "turquoise",
+                    "#E5E5E5",
+                    "#E5E5E5",
+                    "#EEEE00",
+                    "#EEC900",
+                    "darkorange",
+                    "orangered",
+                    "red",
+                    "firebrick",
+                    "mediumvioletred",
+                    "darkviolet",
+                ]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+                cs_1 = plt.pcolormesh(
+                    lon_shift, lat_shift, uh25, transform=transform, cmap=cm, norm=norm
+                )
+                cs_1.cmap.set_under("darkblue")
+                cs_1.cmap.set_over("black")
+                cbar1 = plt.colorbar(
+                    cs_1, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+                )
+                cbar1.set_label(units, fontsize=8)
+                cbar1.ax.tick_params(labelsize=8)
+                ax.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=8,
+                    transform=ax.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT + "/uh25_" + dom + "_f" + fhour + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot Max/Min Hourly 2-5 km UH for: " + dom) % t3)
+    
+            ######################################################
+    
+            t3dom = round(t2 - t1dom, 3)
+            logging.info(("%.3f seconds to plot all variables for forecast hour " + fhour) % t3dom)
+            plt.clf()
+    
+        ######################################################
+    
+        main()
diff --git a/scripts/exregional_plot_allvars_diff.py b/scripts/exregional_plot_allvars_diff.py
new file mode 100755
index 0000000000..e94ec7717f
--- /dev/null
+++ b/scripts/exregional_plot_allvars_diff.py
@@ -0,0 +1,2396 @@
+#!/usr/bin/env python3
+
+################################################################################
+####  Python Script Documentation Block
+#
+# Script name:       	plot_allvars_diff.py
+# Script description:  	Generates difference plots from FV3-LAM post processed
+#                       grib2 output over the CONUS
+#
+# Authors:  Ben Blake		Org: NOAA/NWS/NCEP/EMC		Date: 2020-08-24
+#           David Wright 	Org: University of Michigan
+#
+# Instructions:		Make sure all the necessary modules can be imported.
+#                       The following command line arguments are needed:
+#                       1. Cycle date/time in YYYYMMDDHH format
+#                       2. Starting forecast hour
+#                       3. Ending forecast hour
+#                       4. Forecast hour increment
+#                       5. COMOUT_1: COMOUT directory containing postprocessed data.
+#                       6. COMOUT_2: COMOUT directory for second experiment
+#                       7. CARTOPY_DIR:  Base directory of cartopy shapefiles
+#                          -Shapefiles cannot be directly downloaded to NOAA
+#                            machines from the internet, so shapefiles need to
+#                            be downloaded if geopolitical boundaries are
+#                            desired on the maps.
+#                          -File structure should be:
+#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
+#                          -More information regarding files needed to setup
+#                            display maps in Cartopy, see SRW App Users' Guide
+#                       8. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
+#                          used in forecasts and in constructing the names
+#                          of the post output files.  This must be the same
+#                          for both forecasts.
+#
+#           		To create difference plots for two forecasts on the
+#                       RRFS_CONUS_25km grid for hours 20-24 from the 5/7
+#                       00Z cycle with hourly output:
+#                          python plot_allvars_diff.py 2020050700  20 24 1 \
+#                          /path/to/expt_dir_1 /path/to/expt_dir_2 \
+#                          /path/to/base/cartopy/maps RRFS_CONUS_25km
+#
+#                       The variable domains in this script can be set to either
+#                         'conus' for a CONUS map or 'regional' where the map
+#                         is defined from variables in the grib2 files
+#
+################################################################################
+
+# -------------Import modules --------------------------#
+import pygrib
+import cartopy.crs as ccrs
+from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
+import cartopy.feature as cfeature
+import matplotlib
+
+matplotlib.use("Agg")
+import io
+import matplotlib.pyplot as plt
+import dateutil.relativedelta, dateutil.parser
+from PIL import Image
+from matplotlib.gridspec import GridSpec
+import numpy as np
+import time, os, sys, multiprocessing
+import multiprocessing.pool
+from scipy import ndimage
+import pyproj
+import argparse
+import cartopy
+import logging
+import warnings
+
+# --------------Define some functions ------------------#
+
+
+def ndate(cdate, hours):
+    if not isinstance(cdate, str):
+        if isinstance(cdate, int):
+            cdate = str(cdate)
+        else:
+            logging.error("ERROR: Input cdate must be string or integer!")
+            sys.exit(1)
+    if not isinstance(hours, int):
+        if isinstance(hours, str):
+            hours = int(hours)
+        else:
+            logging.error("ERROR: Input delta must be string or integer!")
+            sys.exit(1)
+
+    indate = cdate.strip()
+    hh = indate[8:10]
+    yyyy = indate[0:4]
+    mm = indate[4:6]
+    dd = indate[6:8]
+    # set date/time field
+    parseme = yyyy + " " + mm + " " + dd + " " + hh
+    datetime_cdate = dateutil.parser.parse(parseme)
+    valid = datetime_cdate + dateutil.relativedelta.relativedelta(hours=+hours)
+    vyyyy = str(valid.year)
+    vm = str(valid.month).zfill(2)
+    vd = str(valid.day).zfill(2)
+    vh = str(valid.hour).zfill(2)
+    return vyyyy + vm + vd + vh
+
+
+def clear_plotables(ax, keep_ax_lst, fig):
+    #### - step to clear off old plottables but leave the map info - ####
+    if len(keep_ax_lst) == 0:
+        logging.info(
+            "clear_plotables WARNING keep_ax_lst has length 0. Clearing ALL plottables including map info!"
+        )
+    cur_ax_children = ax.get_children()[:]
+    if len(cur_ax_children) > 0:
+        for a in cur_ax_children:
+            if a not in keep_ax_lst:
+                # if the artist isn't part of the initial set up, remove it
+                a.remove()
+
+
+def compress_and_save(filename):
+    #### - compress and save the image - ####
+    ram = io.BytesIO()
+    plt.savefig(ram, format="png", bbox_inches="tight", dpi=150)
+    ram.seek(0)
+    im = Image.open(ram)
+    im2 = im.convert("RGB")  # .convert('P', palette=Image.ADAPTIVE)
+    im2.save(filename, format="PNG")
+
+
+def cmap_t2m():
+    # Create colormap for 2-m temperature
+    # Modified version of the ncl_t2m colormap
+
+    #fmt off
+    r = np.array(
+        [255,128,0,  70, 51, 0,  255,0, 0,  51, 255,255,255,255,255,171,128,128,36,162,255]
+    )
+    g = np.array(
+        [0,  0,  0,  70, 102,162,255,92,128,185,255,214,153,102,0,  0,  0,  68, 36,162,255]
+    )
+    b = np.array(
+        [255,128,128,255,255,255,255,0, 0,  102,0,  112,0,  0,  0,  56, 0,  68, 36,162,255]
+    )
+    #fmt on
+
+    xsize = np.arange(np.size(r))
+    r = r / 255.0
+    g = g / 255.0
+    b = b / 255.0
+    red = []
+    green = []
+    blue = []
+    for i in range(len(xsize)):
+        xNorm = float(i) / (float(np.size(r)) - 1.0)
+        red.append([xNorm, r[i], r[i]])
+        green.append([xNorm, g[i], g[i]])
+        blue.append([xNorm, b[i], b[i]])
+    colorDict = {"red": red, "green": green, "blue": blue}
+    cmap_t2m_coltbl = matplotlib.colors.LinearSegmentedColormap(
+        "CMAP_T2M_COLTBL", colorDict
+    )
+    return cmap_t2m_coltbl
+
+
+def cmap_q2m():
+    # Create colormap for 2-m dew point temperature
+    r = np.array([255, 179, 96, 128, 0, 0, 51, 0, 0, 0, 133, 51, 70, 0, 128, 128, 180])
+    g = np.array(
+        [255, 179, 96, 128, 92, 128, 153, 155, 155, 255, 162, 102, 70, 0, 0, 0, 0]
+    )
+    b = np.array(
+        [255, 179, 96, 0, 0, 0, 102, 155, 255, 255, 255, 255, 255, 128, 255, 128, 128]
+    )
+    xsize = np.arange(np.size(r))
+    r = r / 255.0
+    g = g / 255.0
+    b = b / 255.0
+    red = []
+    green = []
+    blue = []
+    for i in range(len(xsize)):
+        xNorm = float(i) / (float(np.size(r)) - 1.0)
+        red.append([xNorm, r[i], r[i]])
+        green.append([xNorm, g[i], g[i]])
+        blue.append([xNorm, b[i], b[i]])
+    colorDict = {"red": red, "green": green, "blue": blue}
+    cmap_q2m_coltbl = matplotlib.colors.LinearSegmentedColormap(
+        "CMAP_Q2M_COLTBL", colorDict
+    )
+    cmap_q2m_coltbl.set_over(color="deeppink")
+    return cmap_q2m_coltbl
+
+
+def rotate_wind(true_lat, lov_lon, earth_lons, uin, vin, proj, inverse=False):
+    #  Rotate winds from LCC relative to earth relative (or vice-versa if inverse==true)
+    #   This routine is vectorized and *should* work on any size 2D vg and ug arrays.
+    #   Program will quit if dimensions are too large.
+    #
+    # Input args:
+    #  true_lat = True latitidue for LCC projection (single value in degrees)
+    #  lov_lon  = The LOV value from grib (e.g. - -95.0) (single value in degrees)
+    #              Grib doc says: "Lov = orientation of the grid; i.e. the east longitude value of
+    #                              the meridian which is parallel to the Y-axis (or columns of the grid)
+    #                              along which latitude increases as the Y-coordinate increases (the
+    #                              orientation longitude may or may not appear on a particular grid).
+    #
+    #  earth_lons = Earth relative longitudes (can be an array, in degrees)
+    #  uin, vin     = Input winds to rotate
+    #
+    # Returns:
+    #  uout, vout = Output, rotated winds
+    # -----------------------------------------------------------------------------------------------------
+
+    # Get size and length of input u winds, if not 2d, raise an error
+    q = np.shape(uin)
+    ndims = len(q)
+    if ndims > 2:
+        # Raise error and quit!
+        raise SystemExit("Input winds for rotation have greater than 2 dimensions!")
+    if lov_lon > 0.0:
+        lov_lon = lov_lon - 360.0
+    dtr = np.pi / 180.0  # Degrees to radians
+
+    if not isinstance(inverse, bool):
+        raise TypeError("**kwarg inverse must be of type bool.")
+
+    # Compute rotation constant which is also
+    # known as the Lambert cone constant.  In the case
+    # of a polar stereographic projection, this is one.
+    # See the following pdf for excellent documentation
+    # http://www.dtcenter.org/met/users/docs/write_ups/velocity.pdf
+    if proj.lower() == "lcc":
+        rotcon_p = np.sin(true_lat * dtr)
+    elif proj.lower() in ["stere", "spstere", "npstere"]:
+        rotcon_p = 1.0
+    else:
+        raise SystemExit(
+            "Unsupported map projection: " + proj.lower() + " for wind rotation."
+        )
+
+    angles = rotcon_p * (earth_lons - lov_lon) * dtr
+    sinx2 = np.sin(angles)
+    cosx2 = np.cos(angles)
+
+    # Steps below are elementwise products, not matrix mutliplies
+    if inverse == False:
+        # Return the earth relative winds
+        uout = cosx2 * uin + sinx2 * vin
+        vout = -sinx2 * uin + cosx2 * vin
+    elif inverse == True:
+        # Return the grid relative winds
+        uout = cosx2 * uin - sinx2 * vin
+        vout = sinx2 * uin + cosx2 * vin
+
+    return uout, vout
+
+
+def setup_logging(debug=False):
+
+    """Calls initialization functions for logging package, and sets the
+    user-defined level for logging in the script."""
+
+    level = logging.INFO
+    if debug:
+        level = logging.DEBUG
+
+    logging.basicConfig(format="%(levelname)s: %(message)s ", level=level)
+    if debug:
+        logging.info("Logging level set to DEBUG")
+
+# -------------Start of script -------------------------#
+if __name__ == "__main__":
+
+    # Define required positional arguments
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--cycle",
+        "-c",
+        help="Cycle date/time in YYYYMMDDHH format.",
+        required=True,
+    )
+    parser.add_argument(
+        "--start",
+        "-s",
+        help="Starting forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--end",
+        "-e",
+        help="Ending forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--inc",
+        "-i",
+        help="Increment forecast hour.",
+        required=True,
+    )
+    parser.add_argument(
+        "--comout-1",
+        help="Path to directory 1 containing post-processed files.",
+        required=True,
+    )
+    parser.add_argument(
+        "--comout-2",
+        help="Path to directory 2 containing post-processed files.",
+        required=True,
+    )
+    parser.add_argument(
+        "--cartopy-dir",
+        help="Path to base directory of cartopy shapefiles.",
+        required=True,
+    )
+    parser.add_argument(
+        "--domain",
+        "-d",
+        help="Name of native domain used in forecast (and in constructing post file names).",
+        required=True,
+    )
+    parser.add_argument(
+        "--debug",
+        action="store_true",
+        help="Print debug messages",
+    )
+    args = parser.parse_args()
+    
+    setup_logging(args.debug)
+
+    # Throw away python warnings (mostly depreciation.)
+    warnings.simplefilter("ignore")
+
+    # Read date/time, forecast hour, and directory paths from command line
+    ymdh = str(args.cycle)
+    ymd = ymdh[0:8]
+    year = int(ymdh[0:4])
+    month = int(ymdh[4:6])
+    day = int(ymdh[6:8])
+    hour = int(ymdh[8:10])
+    cyc = str(hour).zfill(2)
+    logging.info(f"{year}, {month}, {day}, {hour}")
+    
+    # Define the range of forecast hours to create plots for
+    start_fhr = int(args.start)
+    end_fhr = int(args.end)
+    increment_fhr = int(args.inc)
+    if (start_fhr == end_fhr) or (increment_fhr == 0):
+        fhours = [start_fhr]
+    else:
+        num = int(((end_fhr - start_fhr) / increment_fhr) + 1)
+        fhours = np.linspace(start_fhr, end_fhr, num, dtype="int")
+    logging.info(fhours)
+    
+    COMOUT_1 = str(args.comout_1)
+    COMOUT_2 = str(args.comout_2)
+    CARTOPY_DIR = str(args.cartopy_dir)
+    POST_OUTPUT_DOMAIN_NAME = str(args.domain).lower()
+    
+    # Loop over forecast hours
+    for fhr in fhours:
+        fhour = str(fhr).zfill(3)
+        logging.info("Working on forecast hour " + fhour)
+        itime = ymdh
+        vtime = ndate(itime, int(fhr))
+    
+        # Define the location of the input files
+        data1 = pygrib.open(
+            COMOUT_1
+            + "/rrfs.t"
+            + cyc
+            + "z.prslev.f"
+            + fhour
+            + "."
+            + POST_OUTPUT_DOMAIN_NAME
+            + ".grib2"
+        )
+        data2 = pygrib.open(
+            COMOUT_2
+            + "/rrfs.t"
+            + cyc
+            + "z.prslev.f"
+            + fhour
+            + "."
+            + POST_OUTPUT_DOMAIN_NAME
+            + ".grib2"
+        )
+    
+        # Get the lats and lons
+        grids = [data1, data2]
+        lats = []
+        lons = []
+        lats_shift = []
+        lons_shift = []
+    
+        for data in grids:
+            # Unshifted grid for contours and wind barbs
+            lat, lon = data[1].latlons()
+            lats.append(lat)
+            lons.append(lon)
+    
+            # Shift grid for pcolormesh
+            lat1 = data[1]["latitudeOfFirstGridPointInDegrees"]
+            lon1 = data[1]["longitudeOfFirstGridPointInDegrees"]
+            try:
+                nx = data[1]["Nx"]
+                ny = data[1]["Ny"]
+            except:
+                nx = data[1]["Ni"]
+                ny = data[1]["Nj"]
+            dx = data[1]["DxInMetres"]
+            dy = data[1]["DyInMetres"]
+            pj = pyproj.Proj(data[1].projparams)
+            llcrnrx, llcrnry = pj(lon1, lat1)
+            llcrnrx = llcrnrx - (dx / 2.0)
+            llcrnry = llcrnry - (dy / 2.0)
+            x = llcrnrx + dx * np.arange(nx)
+            y = llcrnry + dy * np.arange(ny)
+            x, y = np.meshgrid(x, y)
+            lon, lat = pj(x, y, inverse=True)
+            lats_shift.append(lat)
+            lons_shift.append(lon)
+    
+        # Unshifted lat/lon arrays grabbed directly using latlons() method
+        lat = lats[0]
+        lon = lons[0]
+        lat2 = lats[1]
+        lon2 = lons[1]
+    
+        # Shifted lat/lon arrays for pcolormesh
+        lat_shift = lats_shift[0]
+        lon_shift = lons_shift[0]
+        lat2_shift = lats_shift[1]
+        lon2_shift = lons_shift[1]
+    
+        Lat0 = data1[1]["LaDInDegrees"]
+        Lon0 = data1[1]["LoVInDegrees"]
+        logging.info(Lat0)
+        logging.info(Lon0)
+    
+        # Specify plotting domains
+        # User can add domains here, just need to specify lat/lon information below
+        # (if dom == 'conus' block)
+        domains = ["conus"]  # Other option is 'regional'
+    
+        ###################################################
+        # Read in all variables and calculate differences #
+        ###################################################
+        t1a = time.perf_counter()
+    
+        # Sea level pressure
+        slp_1 = data1.select(name="Pressure reduced to MSL")[0].values * 0.01
+        slpsmooth_1 = ndimage.gaussian_filter(slp_1, 13.78)
+        slp_2 = data2.select(name="Pressure reduced to MSL")[0].values * 0.01
+        slpsmooth_2 = ndimage.gaussian_filter(slp_2, 13.78)
+        slp_diff = slp_2 - slp_1
+    
+        # 2-m temperature
+        tmp2m_1 = data1.select(name="2 metre temperature")[0].values
+        tmp2m_1 = (tmp2m_1 - 273.15) * 1.8 + 32.0
+        tmp2m_2 = data2.select(name="2 metre temperature")[0].values
+        tmp2m_2 = (tmp2m_2 - 273.15) * 1.8 + 32.0
+        tmp2m_diff = tmp2m_2 - tmp2m_1
+    
+        # 2-m dew point temperature
+        dew2m_1 = data1.select(name="2 metre dewpoint temperature")[0].values
+        dew2m_1 = (dew2m_1 - 273.15) * 1.8 + 32.0
+        dew2m_2 = data2.select(name="2 metre dewpoint temperature")[0].values
+        dew2m_2 = (dew2m_2 - 273.15) * 1.8 + 32.0
+        dew2m_diff = dew2m_2 - dew2m_1
+    
+        # 10-m wind speed
+        uwind_1 = data1.select(name="10 metre U wind component")[0].values * 1.94384
+        vwind_1 = data1.select(name="10 metre V wind component")[0].values * 1.94384
+        uwind_2 = data2.select(name="10 metre U wind component")[0].values * 1.94384
+        vwind_2 = data2.select(name="10 metre V wind component")[0].values * 1.94384
+        # Rotate winds from grid relative to Earth relative
+        uwind_1, vwind_1 = rotate_wind(
+            Lat0, Lon0, lon, uwind_1, vwind_1, "lcc", inverse=False
+        )
+        uwind_2, vwind_2 = rotate_wind(
+            Lat0, Lon0, lon2, uwind_2, vwind_2, "lcc", inverse=False
+        )
+        wspd10m_1 = np.sqrt(uwind_1**2 + vwind_1**2)
+        wspd10m_2 = np.sqrt(uwind_2**2 + vwind_2**2)
+        wspd10m_diff = wspd10m_2 - wspd10m_1
+    
+        # Surface-based CAPE
+        cape_1 = data1.select(
+            name="Convective available potential energy", typeOfLevel="surface"
+        )[0].values
+        cape_2 = data2.select(
+            name="Convective available potential energy", typeOfLevel="surface"
+        )[0].values
+        cape_diff = cape_2 - cape_1
+    
+        # Surface-based CIN
+        cin_1 = data1.select(name="Convective inhibition", typeOfLevel="surface")[0].values
+        cin_2 = data2.select(name="Convective inhibition", typeOfLevel="surface")[0].values
+        cin_diff = cin_2 - cin_1
+    
+        # 500 mb height, wind, vorticity
+        try:
+            z500_1 = data1.select(name="Geopotential Height", level=500)[0].values * 0.1
+            z500_1 = ndimage.gaussian_filter(z500_1, 6.89)
+            z500_2 = data2.select(name="Geopotential Height", level=500)[0].values * 0.1
+            z500_2 = ndimage.filters.gaussian_filter(z500_2, 6.89)
+            z500_diff = z500_2 - z500_1
+            vort500_1 = data1.select(name="Absolute vorticity", level=500)[0].values * 100000
+            vort500_1 = ndimage.filters.gaussian_filter(vort500_1, 1.7225)
+            vort500_1[vort500_1 > 1000] = 0  # Mask out undefined values on domain edge
+            vort500_2 = data2.select(name="Absolute vorticity", level=500)[0].values * 100000
+            vort500_2 = ndimage.filters.gaussian_filter(vort500_2, 1.7225)
+            vort500_2[vort500_2 > 1000] = 0  # Mask out undefined values on domain edge
+            u500_1 = data1.select(name="U component of wind", level=500)[0].values * 1.94384
+            u500_2 = data2.select(name="U component of wind", level=500)[0].values * 1.94384
+            v500_1 = data1.select(name="V component of wind", level=500)[0].values * 1.94384
+            v500_2 = data2.select(name="V component of wind", level=500)[0].values * 1.94384
+            # Rotate winds from grid relative to Earth relative
+            u500_1, v500_1 = rotate_wind(Lat0, Lon0, lon, u500_1, v500_1, "lcc", inverse=False)
+            u500_2, v500_2 = rotate_wind(Lat0, Lon0, lon2, u500_2, v500_2, "lcc", inverse=False)
+        except:
+            u500_1 = None
+            u500_2 = None
+    
+        # 250 mb winds
+        u250_1 = data1.select(name="U component of wind", level=250)[0].values * 1.94384
+        u250_2 = data2.select(name="U component of wind", level=250)[0].values * 1.94384
+        v250_1 = data1.select(name="V component of wind", level=250)[0].values * 1.94384
+        v250_2 = data2.select(name="V component of wind", level=250)[0].values * 1.94384
+        # Rotate winds from grid relative to Earth relative
+        u250_1, v250_1 = rotate_wind(Lat0, Lon0, lon, u250_1, v250_1, "lcc", inverse=False)
+        u250_2, v250_2 = rotate_wind(Lat0, Lon0, lon2, u250_2, v250_2, "lcc", inverse=False)
+        wspd250_1 = np.sqrt(u250_1**2 + v250_1**2)
+        wspd250_2 = np.sqrt(u250_2**2 + v250_2**2)
+        wspd250_diff = wspd250_2 - wspd250_1
+    
+        # Total precipitation
+        qpf_1 = (
+            data1.select(name="Total Precipitation", lengthOfTimeRange=fhr)[0].values
+            * 0.0393701
+        )
+        qpf_2 = (
+            data2.select(name="Total Precipitation", lengthOfTimeRange=fhr)[0].values
+            * 0.0393701
+        )
+        qpf_diff = qpf_2 - qpf_1
+    
+        # Composite reflectivity
+        refc_1 = data1.select(name="Maximum/Composite radar reflectivity")[0].values
+        refc_2 = data2.select(name="Maximum/Composite radar reflectivity")[0].values
+    
+        if fhr > 0:
+            # Max/Min Hourly 2-5 km Updraft Helicity
+            maxuh25_1 = data1.select(
+                stepType="max", parameterName="199", topLevel=5000, bottomLevel=2000
+            )[0].values
+            maxuh25_2 = data2.select(
+                stepType="max", parameterName="199", topLevel=5000, bottomLevel=2000
+            )[0].values
+            minuh25_1 = data1.select(
+                stepType="min", parameterName="200", topLevel=5000, bottomLevel=2000
+            )[0].values
+            minuh25_2 = data2.select(
+                stepType="min", parameterName="200", topLevel=5000, bottomLevel=2000
+            )[0].values
+            maxuh25_1[maxuh25_1 < 10] = 0
+            maxuh25_2[maxuh25_2 < 10] = 0
+            minuh25_1[minuh25_1 > -10] = 0
+            minuh25_2[minuh25_2 > -10] = 0
+            uh25_1 = maxuh25_1 + minuh25_1
+            uh25_2 = maxuh25_2 + minuh25_2
+            uh25_diff = uh25_2 - uh25_1
+    
+        t2a = time.perf_counter()
+        t3a = round(t2a - t1a, 3)
+        logging.info(("%.3f seconds to read all messages") % t3a)
+    
+        ########################################
+        #    START PLOTTING FOR EACH DOMAIN    #
+        ########################################
+    
+        def main():
+    
+            # Number of processes must coincide with the number of domains to plot
+            # pool = multiprocessing.Pool(len(domains))
+            # pool.map(plot_all,domains)
+    
+            # To avoid import multiprocessing recursively on MacOS/Windows etc.
+            # Anyway since we only have one domain for SRW application
+            for dom in domains:
+                plot_all(dom)
+    
+        def plot_all(dom):
+    
+            t1dom = time.perf_counter()
+    
+            # Map corners for each domain
+            if dom == "conus":
+                llcrnrlon = -120.5
+                llcrnrlat = 21.0
+                urcrnrlon = -64.5
+                urcrnrlat = 49.0
+                lat_0 = 35.4
+                lon_0 = -97.6
+                extent = [llcrnrlon - 3, urcrnrlon - 6, llcrnrlat - 1, urcrnrlat + 2]
+            elif dom == "regional":
+                llcrnrlon = np.min(lon)
+                llcrnrlat = np.min(lat)
+                urcrnrlon = np.max(lon)
+                urcrnrlat = np.max(lat)
+                lat_0 = Lat0
+                lon_0 = Lon0
+                extent = [llcrnrlon, urcrnrlon, llcrnrlat - 1, urcrnrlat]
+    
+            # create figure and axes instances
+            fig = plt.figure(figsize=(10, 10))
+            gs = GridSpec(9, 9, wspace=0.0, hspace=0.0)
+    
+            # Define where Cartopy Maps are located
+            cartopy.config["data_dir"] = CARTOPY_DIR
+    
+            back_res = "50m"
+            back_img = "on"
+    
+            # set up the map background with cartopy
+            myproj = ccrs.LambertConformal(
+                central_longitude=lon_0,
+                central_latitude=lat_0,
+                false_easting=0.0,
+                false_northing=0.0,
+                secant_latitudes=None,
+                standard_parallels=None,
+                globe=None,
+            )
+            ax1 = fig.add_subplot(gs[0:4, 0:4], projection=myproj)
+            ax2 = fig.add_subplot(gs[0:4, 5:], projection=myproj)
+            ax3 = fig.add_subplot(gs[5:, 1:8], projection=myproj)
+            ax1.set_extent(extent)
+            ax2.set_extent(extent)
+            ax3.set_extent(extent)
+    
+            fline_wd = 0.5  # line width
+            falpha = 0.3  # transparency
+    
+            # natural_earth
+            #    land=cfeature.NaturalEarthFeature('physical','land',back_res,
+            #                      edgecolor='face',facecolor=cfeature.COLORS['land'],
+            #                      alpha=falpha)
+            lakes = cfeature.NaturalEarthFeature(
+                "physical",
+                "lakes",
+                back_res,
+                edgecolor="blue",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+            coastline = cfeature.NaturalEarthFeature(
+                "physical",
+                "coastline",
+                back_res,
+                edgecolor="blue",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+            states = cfeature.NaturalEarthFeature(
+                "cultural",
+                "admin_1_states_provinces",
+                back_res,
+                edgecolor="black",
+                facecolor="none",
+                linewidth=fline_wd,
+                linestyle=":",
+                alpha=falpha,
+            )
+            borders = cfeature.NaturalEarthFeature(
+                "cultural",
+                "admin_0_countries",
+                back_res,
+                edgecolor="red",
+                facecolor="none",
+                linewidth=fline_wd,
+                alpha=falpha,
+            )
+    
+            # All lat lons are earth relative, so setup the associated projection correct for that data
+            transform = ccrs.PlateCarree()
+    
+            # high-resolution background images
+            if back_img == "on":
+                img = plt.imread(CARTOPY_DIR + "/raster_files/NE1_50M_SR_W.tif")
+                ax1.imshow(img, origin="upper", transform=transform)
+                ax2.imshow(img, origin="upper", transform=transform)
+                ax3.imshow(img, origin="upper", transform=transform)
+    
+            #  ax.add_feature(land)
+            ax1.add_feature(lakes)
+            ax1.add_feature(states)
+            ax1.add_feature(borders)
+            ax1.add_feature(coastline)
+            ax2.add_feature(lakes)
+            ax2.add_feature(states)
+            ax2.add_feature(borders)
+            ax2.add_feature(coastline)
+            ax3.add_feature(lakes)
+            ax3.add_feature(states)
+            ax3.add_feature(borders)
+            ax3.add_feature(coastline)
+    
+            # Map/figure has been set up here, save axes instances for use again later
+            keep_ax_lst_1 = ax1.get_children()[:]
+            keep_ax_lst_2 = ax2.get_children()[:]
+            keep_ax_lst_3 = ax3.get_children()[:]
+    
+            # colors for difference plots, only need to define once
+            diffcolors = [
+                "blue",
+                "#1874CD",
+                "dodgerblue",
+                "deepskyblue",
+                "turquoise",
+                "white",
+                "white",
+                "#EEEE00",
+                "#EEC900",
+                "darkorange",
+                "orangered",
+                "red",
+            ]
+    
+            ################################
+            # Plot SLP
+            ################################
+            t1 = time.perf_counter()
+            logging.info(("Working on slp for " + dom))
+    
+            units = "mb"
+            clevs = [
+                976,
+                980,
+                984,
+                988,
+                992,
+                996,
+                1000,
+                1004,
+                1008,
+                1012,
+                1016,
+                1020,
+                1024,
+                1028,
+                1032,
+                1036,
+                1040,
+                1044,
+                1048,
+                1052,
+            ]
+            clevsdiff = [-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12]
+            cm = plt.cm.Spectral_r
+            cmdiff = matplotlib.colors.ListedColormap(diffcolors)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs1_a = ax1.pcolormesh(
+                lon_shift, lat_shift, slp_1, transform=transform, cmap=cm, norm=norm
+            )
+            cbar1 = plt.colorbar(
+                cs1_a, ax=ax1, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=5)
+            cs1_b = ax1.contour(
+                lon_shift,
+                lat_shift,
+                slpsmooth_1,
+                np.arange(940, 1060, 4),
+                colors="black",
+                linewidths=1.25,
+                transform=transform,
+            )
+            plt.clabel(cs1_b, np.arange(940, 1060, 4), inline=1, fmt="%d", fontsize=6)
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM SLP ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs2_a = ax2.pcolormesh(
+                lon2_shift, lat2_shift, slp_2, transform=transform, cmap=cm, norm=norm
+            )
+            cbar2 = plt.colorbar(
+                cs2_a, ax=ax2, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=5)
+            cs2_b = ax2.contour(
+                lon2_shift,
+                lat2_shift,
+                slpsmooth_2,
+                np.arange(940, 1060, 4),
+                colors="black",
+                linewidths=1.25,
+                transform=transform,
+            )
+            plt.clabel(cs2_b, np.arange(940, 1060, 4), inline=1, fmt="%d", fontsize=6)
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 SLP ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=8,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                slp_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=5)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM SLP ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1 + "/slp_diff_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot slp for: " + dom) % t3)
+    
+            #################################
+            # Plot 2-m T
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on t2m for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "\xb0" "F"
+            clevs = np.linspace(-16, 134, 51)
+            clevsdiff = [-6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6]
+            #  cm = plt.cm.Spectral_r
+            cm = cmap_t2m()
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift, lat_shift, tmp2m_1, transform=transform, cmap=cm, norm=norm
+            )
+            cs_1.cmap.set_under("white")
+            cs_1.cmap.set_over("white")
+            cbar1 = plt.colorbar(
+                cs_1,
+                ax=ax1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=[-16, -4, 8, 20, 32, 44, 56, 68, 80, 92, 104, 116, 128],
+                extend="both",
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=5)
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM 2-m Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift, lat2_shift, tmp2m_2, transform=transform, cmap=cm, norm=norm
+            )
+            cs_2.cmap.set_under("white")
+            cs_2.cmap.set_over("white")
+            cbar2 = plt.colorbar(
+                cs_2,
+                ax=ax2,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=[-16, -4, 8, 20, 32, 44, 56, 68, 80, 92, 104, 116, 128],
+                extend="both",
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=5)
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 2-m Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                tmp2m_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=6)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM 2-m Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1 + "/2mt_diff_" + dom + "_f" + fhour + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 2mt for: " + dom) % t3)
+    
+            #################################
+            # Plot 2-m Dew Point
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 2mdew for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "\xb0" "F"
+            clevs = np.linspace(-5, 80, 35)
+            clevsdiff = [-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12]
+            cm = cmap_q2m()
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift, lat_shift, dew2m_1, transform=transform, cmap=cm, norm=norm
+            )
+            cbar1 = plt.colorbar(
+                cs_1, ax=ax1, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=6)
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM 2-m Dew Point Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift, lat2_shift, dew2m_2, transform=transform, cmap=cm, norm=norm
+            )
+            cbar2 = plt.colorbar(
+                cs_2, ax=ax2, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=6)
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 2-m Dew Point Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                dew2m_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=6)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM 2-m Dew Point Temperature ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1
+                + "/2mdew_diff_"
+                + dom
+                + "_f"
+                + fhour
+                + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 2mdew for: " + dom) % t3)
+    
+            #################################
+            # Plot 10-m WSPD
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 10mwspd for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "kts"
+            # Places a wind barb every ~180 km, optimized for CONUS domain
+            skip = round(177.28 * (dx / 1000.0) ** -0.97)
+            logging.info("skipping every " + str(skip) + " grid points to plot")
+            barblength = 4
+    
+            clevs = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60]
+            clevsdiff = [-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12]
+            colorlist = [
+                "turquoise",
+                "dodgerblue",
+                "blue",
+                "#FFF68F",
+                "#E3CF57",
+                "peru",
+                "brown",
+                "crimson",
+                "red",
+                "fuchsia",
+                "DarkViolet",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift,
+                lat_shift,
+                wspd10m_1,
+                transform=transform,
+                cmap=cm,
+                vmin=5,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                ax=ax1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=6)
+            ax1.barbs(
+                lon_shift[::skip, ::skip],
+                lat_shift[::skip, ::skip],
+                uwind_1[::skip, ::skip],
+                vwind_1[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM 10-m Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                wspd10m_2,
+                transform=transform,
+                cmap=cm,
+                vmin=5,
+                norm=norm,
+            )
+            cs_2.cmap.set_under("white", alpha=0.0)
+            cs_2.cmap.set_over("black")
+            cbar2 = plt.colorbar(
+                cs_2,
+                ax=ax2,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=6)
+            ax2.barbs(
+                lon2_shift[::skip, ::skip],
+                lat2_shift[::skip, ::skip],
+                uwind_2[::skip, ::skip],
+                vwind_2[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 10-m Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                wspd10m_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=6)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM 10-m Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1
+                + "/10mwind_diff_"
+                + dom
+                + "_f"
+                + fhour
+                + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 10mwspd for: " + dom) % t3)
+    
+            #################################
+            # Plot Surface-Based CAPE/CIN
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on surface-based CAPE/CIN for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "J/kg"
+            clevs = [100, 250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000]
+            clevs2 = [-2000, -500, -250, -100, -25]
+            clevsdiff = [
+                -2000,
+                -1500,
+                -1000,
+                -500,
+                -250,
+                -100,
+                0,
+                100,
+                250,
+                500,
+                1000,
+                1500,
+                2000,
+            ]
+            colorlist = [
+                "blue",
+                "dodgerblue",
+                "cyan",
+                "mediumspringgreen",
+                "#FAFAD2",
+                "#EEEE00",
+                "#EEC900",
+                "darkorange",
+                "crimson",
+                "darkred",
+                "darkviolet",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift,
+                lat_shift,
+                cape_1,
+                transform=transform,
+                cmap=cm,
+                vmin=100,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                ax=ax1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=4)
+            cs_1b = ax1.contourf(
+                lon_shift,
+                lat_shift,
+                cin_1,
+                clevs2,
+                colors="none",
+                hatches=["**", "++", "////", ".."],
+                transform=transform,
+            )
+            ax1.text(
+                0.5,
+                1.05,
+                "FV3-LAM Surface-Based CAPE (shaded) and CIN (hatched) ("
+                + units
+                + ") \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                cape_2,
+                transform=transform,
+                cmap=cm,
+                vmin=100,
+                norm=norm,
+            )
+            cs_2.cmap.set_under("white", alpha=0.0)
+            cs_2.cmap.set_over("black")
+            cbar2 = plt.colorbar(
+                cs_2,
+                ax=ax2,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=4)
+            cs_2b = ax2.contourf(
+                lon2_shift,
+                lat2_shift,
+                cin_2,
+                clevs2,
+                colors="none",
+                hatches=["**", "++", "////", ".."],
+                transform=transform,
+            )
+            ax2.text(
+                0.5,
+                1.05,
+                "FV3-LAM-2 Surface-Based CAPE (shaded) and CIN (hatched) ("
+                + units
+                + ") \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                cape_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=6)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM Surface-Based CAPE ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1
+                + "/sfcape_diff_"
+                + dom
+                + "_f"
+                + fhour
+                + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot surface-based CAPE/CIN for: " + dom) % t3)
+    
+            #################################
+            # Plot 500 mb HGT/WIND/VORT
+            #################################
+            if (u500_1 is not None) and (u500_2 is not None):
+                t1 = time.perf_counter()
+                logging.info(("Working on 500 mb Hgt/Wind/Vort for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                cbar2.remove()
+                cbar3.remove()
+                clear_plotables(ax1, keep_ax_lst_1, fig)
+                clear_plotables(ax2, keep_ax_lst_2, fig)
+                clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+                units = "x10${^5}$ s${^{-1}}$"
+                skip = round(177.28 * (dx / 1000.0) ** -0.97)
+                logging.info("skipping every " + str(skip) + " grid points to plot")
+                barblength = 4
+    
+                vortlevs = [16, 20, 24, 28, 32, 36, 40]
+                clevsdiff = [-6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6]
+                colorlist = ["yellow", "gold", "goldenrod", "orange", "orangered", "red"]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                norm = matplotlib.colors.BoundaryNorm(vortlevs, cm.N)
+                normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+                cs1_a = ax1.pcolormesh(
+                    lon_shift, lat_shift, vort500_1, transform=transform, cmap=cm, norm=norm
+                )
+                cs1_a.cmap.set_under("white")
+                cs1_a.cmap.set_over("darkred")
+                cbar1 = plt.colorbar(
+                    cs1_a,
+                    ax=ax1,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    ticks=vortlevs,
+                    extend="both",
+                )
+                cbar1.set_label(units, fontsize=6)
+                cbar1.ax.tick_params(labelsize=6)
+                ax1.barbs(
+                    lon_shift[::skip, ::skip],
+                    lat_shift[::skip, ::skip],
+                    u500_1[::skip, ::skip],
+                    v500_1[::skip, ::skip],
+                    length=barblength,
+                    linewidth=0.5,
+                    color="steelblue",
+                    transform=transform,
+                )
+                cs1_b = ax1.contour(
+                    lon_shift,
+                    lat_shift,
+                    z500_1,
+                    np.arange(486, 600, 6),
+                    colors="black",
+                    linewidths=1,
+                    transform=transform,
+                )
+                plt.clabel(
+                    cs1_b, np.arange(486, 600, 6), inline_spacing=1, fmt="%d", fontsize=8
+                )
+                ax1.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM 500 mb Heights (dam), Winds (kts), and $\zeta$ ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=8,
+                    transform=ax1.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs2_a = ax2.pcolormesh(
+                    lon2_shift, lat2_shift, vort500_2, transform=transform, cmap=cm, norm=norm
+                )
+                cs2_a.cmap.set_under("white")
+                cs2_a.cmap.set_over("darkred")
+                cbar2 = plt.colorbar(
+                    cs2_a,
+                    ax=ax2,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    ticks=vortlevs,
+                    extend="both",
+                )
+                cbar2.set_label(units, fontsize=6)
+                cbar2.ax.tick_params(labelsize=6)
+                ax2.barbs(
+                    lon2_shift[::skip, ::skip],
+                    lat2_shift[::skip, ::skip],
+                    u500_2[::skip, ::skip],
+                    v500_2[::skip, ::skip],
+                    length=barblength,
+                    linewidth=0.5,
+                    color="steelblue",
+                    transform=transform,
+                )
+                cs2_b = ax2.contour(
+                    lon2_shift,
+                    lat2_shift,
+                    z500_2,
+                    np.arange(486, 600, 6),
+                    colors="black",
+                    linewidths=1,
+                    transform=transform,
+                )
+                plt.clabel(
+                    cs2_b, np.arange(486, 600, 6), inline_spacing=1, fmt="%d", fontsize=8
+                )
+                ax2.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 500 mb Heights (dam), Winds (kts), and $\zeta$ ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=8,
+                    transform=ax2.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs = ax3.pcolormesh(
+                    lon2_shift,
+                    lat2_shift,
+                    z500_diff,
+                    transform=transform,
+                    cmap=cmdiff,
+                    norm=normdiff,
+                )
+                cs.cmap.set_under("darkblue")
+                cs.cmap.set_over("darkred")
+                cbar3 = plt.colorbar(
+                    cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+                )
+                cbar3.set_label(units, fontsize=6)
+                cbar3.ax.tick_params(labelsize=6)
+                ax3.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 - FV3-LAM 500-mb Heights (dam) \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax3.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT_1 + "/500_diff_" + dom + "_f" + fhour + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot 500 mb Hgt/Wind/Vort for: " + dom) % t3)
+    
+            #################################
+            # Plot 250 mb WIND
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on 250 mb WIND for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "kts"
+            skip = round(177.28 * (dx / 1000.0) ** -0.97)
+            logging.info("skipping every " + str(skip) + " grid points to plot")
+            barblength = 4
+    
+            clevs = [50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150]
+            clevsdiff = [-30, -25, -20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30]
+            colorlist = [
+                "turquoise",
+                "deepskyblue",
+                "dodgerblue",
+                "#1874CD",
+                "blue",
+                "beige",
+                "khaki",
+                "peru",
+                "brown",
+                "crimson",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+            normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift,
+                lat_shift,
+                wspd250_1,
+                transform=transform,
+                cmap=cm,
+                vmin=50,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("red")
+            cbar1 = plt.colorbar(
+                cs_1,
+                ax=ax1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=6)
+            ax1.barbs(
+                lon_shift[::skip, ::skip],
+                lat_shift[::skip, ::skip],
+                u250_1[::skip, ::skip],
+                v250_1[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM 250 mb Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                wspd250_2,
+                transform=transform,
+                cmap=cm,
+                vmin=50,
+                norm=norm,
+            )
+            cs_2.cmap.set_under("white", alpha=0.0)
+            cs_2.cmap.set_over("red")
+            cbar2 = plt.colorbar(
+                cs_2,
+                ax=ax2,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=6)
+            ax2.barbs(
+                lon2_shift[::skip, ::skip],
+                lat2_shift[::skip, ::skip],
+                u250_2[::skip, ::skip],
+                v250_2[::skip, ::skip],
+                length=barblength,
+                linewidth=0.5,
+                color="black",
+                transform=transform,
+            )
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 250 mb Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs = ax3.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                wspd250_diff,
+                transform=transform,
+                cmap=cmdiff,
+                norm=normdiff,
+            )
+            cs.cmap.set_under("darkblue")
+            cs.cmap.set_over("darkred")
+            cbar3 = plt.colorbar(
+                cs, ax=ax3, orientation="horizontal", pad=0.05, shrink=0.6, extend="both"
+            )
+            cbar3.set_label(units, fontsize=6)
+            cbar3.ax.tick_params(labelsize=6)
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 - FV3-LAM 250-mb Winds ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1
+                + "/250wind_diff_"
+                + dom
+                + "_f"
+                + fhour
+                + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot 250 mb WIND for: " + dom) % t3)
+    
+            #################################
+            # Plot Total QPF
+            #################################
+            if fhr > 0:  # Do not make total QPF plot for forecast hour 0
+                t1 = time.perf_counter()
+                logging.info(("Working on total qpf for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                cbar2.remove()
+                cbar3.remove()
+                clear_plotables(ax1, keep_ax_lst_1, fig)
+                clear_plotables(ax2, keep_ax_lst_2, fig)
+                clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+                units = "in"
+                clevs = [
+                    0.01,
+                    0.1,
+                    0.25,
+                    0.5,
+                    0.75,
+                    1,
+                    1.25,
+                    1.5,
+                    1.75,
+                    2,
+                    2.5,
+                    3,
+                    4,
+                    5,
+                    7,
+                    10,
+                    15,
+                    20,
+                ]
+                clevsdiff = [-3, -2.5, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3]
+                colorlist = [
+                    "chartreuse",
+                    "limegreen",
+                    "green",
+                    "blue",
+                    "dodgerblue",
+                    "deepskyblue",
+                    "cyan",
+                    "mediumpurple",
+                    "mediumorchid",
+                    "darkmagenta",
+                    "darkred",
+                    "crimson",
+                    "orangered",
+                    "darkorange",
+                    "goldenrod",
+                    "gold",
+                    "yellow",
+                ]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+                normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+                cs_1 = ax1.pcolormesh(
+                    lon_shift,
+                    lat_shift,
+                    qpf_1,
+                    transform=transform,
+                    cmap=cm,
+                    vmin=0.01,
+                    norm=norm,
+                )
+                cs_1.cmap.set_under("white", alpha=0.0)
+                cs_1.cmap.set_over("pink")
+                cbar1 = plt.colorbar(
+                    cs_1,
+                    ax=ax1,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="max",
+                )
+                cbar1.set_label(units, fontsize=6)
+                #cbar1.ax.set_xticklabels([0.1, 0.5, 1, 1.5, 2, 3, 5, 10, 20])
+                cbar1.ax.set_xticklabels([0.1, 0.3, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 7.5, 10, 15, 17.5, 20])
+                cbar1.ax.tick_params(labelsize=6)
+                ax1.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM "
+                    + fhour
+                    + "-hr Accumulated Precipitation ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax1.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs_2 = ax2.pcolormesh(
+                    lon2_shift,
+                    lat2_shift,
+                    qpf_2,
+                    transform=transform,
+                    cmap=cm,
+                    vmin=0.01,
+                    norm=norm,
+                )
+                cs_2.cmap.set_under("white", alpha=0.0)
+                cs_2.cmap.set_over("pink")
+                cbar2 = plt.colorbar(
+                    cs_2,
+                    ax=ax2,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="max",
+                )
+                cbar2.set_label(units, fontsize=6)
+                #cbar2.ax.set_xticklabels([0.1, 0.5, 1, 1.5, 2, 3, 5, 10, 20])
+                cbar2.ax.set_xticklabels([0.1, 0.3, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 4, 5, 7.5, 10, 15, 17.5, 20])
+                cbar2.ax.tick_params(labelsize=6)
+                ax2.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 "
+                    + fhour
+                    + "-hr Accumulated Precipitation ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax2.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs = ax3.pcolormesh(
+                    lon2_shift,
+                    lat2_shift,
+                    qpf_diff,
+                    transform=transform,
+                    cmap=cmdiff,
+                    norm=normdiff,
+                )
+                cs.cmap.set_under("darkblue")
+                cs.cmap.set_over("darkred")
+                cbar3 = plt.colorbar(
+                    cs,
+                    ax=ax3,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="both",
+                )
+                cbar3.set_label(units, fontsize=6)
+                cbar3.ax.tick_params(labelsize=6)
+                ax3.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 - FV3-LAM "
+                    + fhour
+                    + "-hr Accumulated Precipitation ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax3.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT_1
+                    + "/qpf_diff_"
+                    + dom
+                    + "_f"
+                    + fhour
+                    + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot total qpf for: " + dom) % t3)
+    
+                #################################
+                # Plot Max/Min Hourly 2-5 km UH
+                #################################
+                # Do not make max/min hourly 2-5 km UH plot for forecast hour 0
+                t1 = time.perf_counter()
+                logging.info(("Working on Max/Min Hourly 2-5 km UH for " + dom))
+    
+                # Clear off old plottables but keep all the map info
+                cbar1.remove()
+                cbar2.remove()
+                cbar3.remove()
+                clear_plotables(ax1, keep_ax_lst_1, fig)
+                clear_plotables(ax2, keep_ax_lst_2, fig)
+                clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+                units = "m${^2}$ s$^{-2}$"
+                clevs = [
+                    -150,
+                    -100,
+                    -75,
+                    -50,
+                    -25,
+                    -10,
+                    0,
+                    10,
+                    25,
+                    50,
+                    75,
+                    100,
+                    150,
+                    200,
+                    250,
+                    300,
+                ]
+                #   alternative colormap for just max UH if you don't want to plot the min UH too
+                #   colorlist = ['white','skyblue','mediumblue','green','orchid','firebrick','#EEC900','DarkViolet']
+                colorlist = [
+                    "blue",
+                    "#1874CD",
+                    "dodgerblue",
+                    "deepskyblue",
+                    "turquoise",
+                    "#E5E5E5",
+                    "#E5E5E5",
+                    "#EEEE00",
+                    "#EEC900",
+                    "darkorange",
+                    "orangered",
+                    "red",
+                    "firebrick",
+                    "mediumvioletred",
+                    "darkviolet",
+                ]
+                cm = matplotlib.colors.ListedColormap(colorlist)
+                cmdiff = matplotlib.colors.ListedColormap(diffcolors)
+                norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+                normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
+    
+                cs_1 = ax1.pcolormesh(
+                    lon_shift, lat_shift, uh25_1, transform=transform, cmap=cm, norm=norm
+                )
+                cs_1.cmap.set_under("darkblue")
+                cs_1.cmap.set_over("black")
+                cbar1 = plt.colorbar(
+                    cs_1,
+                    ax=ax1,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="both",
+                )
+                cbar1.set_label(units, fontsize=6)
+                cbar1.ax.tick_params(labelsize=6)
+                ax1.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax1.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs_2 = ax2.pcolormesh(
+                    lon2_shift, lat2_shift, uh25_2, transform=transform, cmap=cm, norm=norm
+                )
+                cs_2.cmap.set_under("darkblue")
+                cs_2.cmap.set_over("black")
+                cbar2 = plt.colorbar(
+                    cs_2,
+                    ax=ax2,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="both",
+                )
+                cbar2.set_label(units, fontsize=6)
+                cbar2.ax.tick_params(labelsize=6)
+                ax2.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 1-h Max/Min 2-5 km Updraft Helicity ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax2.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                cs = ax3.pcolormesh(
+                    lon2_shift,
+                    lat2_shift,
+                    uh25_diff,
+                    transform=transform,
+                    cmap=cmdiff,
+                    norm=normdiff,
+                )
+                cs.cmap.set_under("darkblue")
+                cs.cmap.set_over("darkred")
+                cbar3 = plt.colorbar(
+                    cs,
+                    ax=ax3,
+                    orientation="horizontal",
+                    pad=0.05,
+                    shrink=0.6,
+                    extend="both",
+                )
+                cbar3.set_label(units, fontsize=6)
+                cbar3.ax.tick_params(labelsize=6)
+                ax3.text(
+                    0.5,
+                    1.03,
+                    "FV3-LAM-2 - FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ("
+                    + units
+                    + ") \n initialized: "
+                    + itime
+                    + " valid: "
+                    + vtime
+                    + " (f"
+                    + fhour
+                    + ")",
+                    horizontalalignment="center",
+                    fontsize=6,
+                    transform=ax3.transAxes,
+                    bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+                )
+    
+                compress_and_save(
+                    COMOUT_1
+                    + "/uh25_diff_"
+                    + dom
+                    + "_f"
+                    + fhour
+                    + ".png"
+                )
+                t2 = time.perf_counter()
+                t3 = round(t2 - t1, 3)
+                logging.info(("%.3f seconds to plot Max/Min Hourly 2-5 km UH for: " + dom) % t3)
+    
+            #################################
+            # Plot composite reflectivity
+            #################################
+            t1 = time.perf_counter()
+            logging.info(("Working on composite reflectivity for " + dom))
+    
+            # Clear off old plottables but keep all the map info
+            cbar1.remove()
+            cbar2.remove()
+            cbar3.remove()
+            clear_plotables(ax1, keep_ax_lst_1, fig)
+            clear_plotables(ax2, keep_ax_lst_2, fig)
+            clear_plotables(ax3, keep_ax_lst_3, fig)
+    
+            units = "dBZ"
+            clevs = np.linspace(5, 70, 14)
+            clevsdiff = [20, 1000]
+            colorlist = [
+                "turquoise",
+                "dodgerblue",
+                "mediumblue",
+                "lime",
+                "limegreen",
+                "green",
+                "#EEEE00",
+                "#EEC900",
+                "darkorange",
+                "red",
+                "firebrick",
+                "darkred",
+                "fuchsia",
+            ]
+            cm = matplotlib.colors.ListedColormap(colorlist)
+            norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
+    
+            cs_1 = ax1.pcolormesh(
+                lon_shift,
+                lat_shift,
+                refc_1,
+                transform=transform,
+                cmap=cm,
+                vmin=5,
+                norm=norm,
+            )
+            cs_1.cmap.set_under("white", alpha=0.0)
+            cs_1.cmap.set_over("black")
+            cbar1 = plt.colorbar(
+                cs_1,
+                ax=ax1,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar1.set_label(units, fontsize=6)
+            cbar1.ax.tick_params(labelsize=6)
+            ax1.text(
+                0.5,
+                1.03,
+                "FV3-LAM Composite Reflectivity ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax1.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            cs_2 = ax2.pcolormesh(
+                lon2_shift,
+                lat2_shift,
+                refc_2,
+                transform=transform,
+                cmap=cm,
+                vmin=5,
+                norm=norm,
+            )
+            cs_2.cmap.set_under("white", alpha=0.0)
+            cs_2.cmap.set_over("black")
+            cbar2 = plt.colorbar(
+                cs_2,
+                ax=ax2,
+                orientation="horizontal",
+                pad=0.05,
+                shrink=0.6,
+                ticks=clevs,
+                extend="max",
+            )
+            cbar2.set_label(units, fontsize=6)
+            cbar2.ax.tick_params(labelsize=6)
+            ax2.text(
+                0.5,
+                1.03,
+                "FV3-LAM-2 Composite Reflectivity ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax2.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            csdiff = ax3.contourf(
+                lon_shift, lat_shift, refc_1, clevsdiff, colors="red", transform=transform
+            )
+            csdiff2 = ax3.contourf(
+                lon2_shift,
+                lat2_shift,
+                refc_2,
+                clevsdiff,
+                colors="dodgerblue",
+                transform=transform,
+            )
+            ax3.text(
+                0.5,
+                1.03,
+                "FV3-LAM (red) and FV3-LAM-2 (blue) Composite Reflectivity > 20 ("
+                + units
+                + ") \n initialized: "
+                + itime
+                + " valid: "
+                + vtime
+                + " (f"
+                + fhour
+                + ")",
+                horizontalalignment="center",
+                fontsize=6,
+                transform=ax3.transAxes,
+                bbox=dict(facecolor="white", alpha=0.85, boxstyle="square,pad=0.2"),
+            )
+    
+            compress_and_save(
+                COMOUT_1
+                + "/refc_diff_"
+                + dom
+                + "_f"
+                + fhour
+                + ".png"
+            )
+            t2 = time.perf_counter()
+            t3 = round(t2 - t1, 3)
+            logging.info(("%.3f seconds to plot composite reflectivity for: " + dom) % t3)
+    
+            ######################################################
+    
+            t3dom = round(t2 - t1dom, 3)
+            logging.info(("%.3f seconds to plot all variables for: " + dom) % t3dom)
+            plt.clf()
+    
+        ######################################################
+    
+        main()
diff --git a/tests/WE2E/test_configs/grids_extrn_mdls_suites_community/config.grid_RRFS_CONUS_25km_ics_FV3GFS_lbcs_FV3GFS_suite_GFS_2017_gfdlmp_regional_plot.yaml b/tests/WE2E/test_configs/grids_extrn_mdls_suites_community/config.grid_RRFS_CONUS_25km_ics_FV3GFS_lbcs_FV3GFS_suite_GFS_2017_gfdlmp_regional_plot.yaml
new file mode 100644
index 0000000000..1c837100e5
--- /dev/null
+++ b/tests/WE2E/test_configs/grids_extrn_mdls_suites_community/config.grid_RRFS_CONUS_25km_ics_FV3GFS_lbcs_FV3GFS_suite_GFS_2017_gfdlmp_regional_plot.yaml
@@ -0,0 +1,23 @@
+metadata:
+  description: |-
+    This test is to ensure that the workflow running in community mode can produce
+    graphics on the RRFS_CONUS_25km grid and GFS_2017_gfdlmp_regional
+    physics suite with ICs and LBCs derived from the FV3GFS.
+user:
+  RUN_ENVIR: community
+workflow:
+  CCPP_PHYS_SUITE: FV3_GFS_2017_gfdlmp_regional
+  PREDEF_GRID_NAME: RRFS_CONUS_25km
+  DATE_FIRST_CYCL: '2019070100'
+  DATE_LAST_CYCL: '2019070100'
+  FCST_LEN_HRS: 6
+  PREEXISTING_DIR_METHOD: rename
+workflow_switches:
+  RUN_TASK_PLOT_ALLVARS: true
+task_get_extrn_ics:
+  EXTRN_MDL_NAME_ICS: FV3GFS
+  USE_USER_STAGED_EXTRN_FILES: true
+task_get_extrn_lbcs:
+  EXTRN_MDL_NAME_LBCS: FV3GFS
+  LBC_SPEC_INTVL_HRS: 6
+  USE_USER_STAGED_EXTRN_FILES: true
diff --git a/ush/Python/environment.yml b/ush/Python/environment.yml
deleted file mode 100644
index 8e08987ce0..0000000000
--- a/ush/Python/environment.yml
+++ /dev/null
@@ -1,15 +0,0 @@
-name: pygraf
-channels:
-  - conda-forge
-  - defaults
-dependencies:
-  - python=3.7*
-  - basemap=1.2*
-  - basemap-data-hires=1.2.1
-  - pynio=1.5.5
-  - matplotlib=3.2*
-  - metpy=0.12.1
-  - pylint=2.4*
-  - pytest=6.1*
-  - pyyaml=5.3.1
-  - xarray=0.15.1
diff --git a/ush/Python/plot_allvars.py b/ush/Python/plot_allvars.py
deleted file mode 100755
index 87a38ab0e1..0000000000
--- a/ush/Python/plot_allvars.py
+++ /dev/null
@@ -1,811 +0,0 @@
-################################################################################
-####  Python Script Documentation Block
-#
-# Script name:       	plot_allvars.py
-# Script description:  	Generates plots from FV3-LAM post processed grib2 output
-#			over the CONUS
-#
-# Authors:  Ben Blake		Org: NOAA/NWS/NCEP/EMC		Date: 2020-05-07
-#           David Wright 	Org: University of Michigan
-#
-# Instructions:		Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour
-#                       3. Ending forecast hour
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR: Experiment directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR/YYYYMMDDHH/postprd/
-#                       6. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -Shapefiles cannot be directly downloaded to NOAA
-#                            machines from the internet, so shapefiles need to
-#                            be downloaded if geopolitical boundaries are
-#                            desired on the maps.
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                          -More information regarding files needed to setup
-#                            display maps in Cartopy, see SRW App Users' Guide
-#                       7. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecast and in constructing the names 
-#                          of the post output files.
-#
-#           		To create plots for a forecast on the RRFS_CONUS_25km
-#                       grid for hours 20-24 from the 5/7 00Z cycle with 
-#                       hourly output:
-#                          python plot_allvars.py 2020050700 20 24 1 \
-#                          /path/to/expt_dirs/experiment/name \
-#                          /path/to/base/cartopy/maps RRFS_CONUS_25km
-#
-#                       The variable domains in this script can be set to either
-#                         'conus' for a CONUS map or 'regional' where the map
-#                         is defined from variables in the grib2 files
-#
-################################################################################
-
-#-------------Import modules --------------------------#
-import pygrib
-import cartopy.crs as ccrs
-from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
-import cartopy.feature as cfeature
-import matplotlib
-matplotlib.use('Agg')
-import io
-import matplotlib.pyplot as plt
-import dateutil.relativedelta, dateutil.parser
-from PIL import Image
-from matplotlib.gridspec import GridSpec
-import numpy as np
-import time,os,sys,multiprocessing
-import multiprocessing.pool
-from scipy import ndimage
-import pyproj
-import argparse
-import cartopy
-
-#--------------Define some functions ------------------#
-
-def ndate(cdate,hours):
-   if not isinstance(cdate, str):
-     if isinstance(cdate, int):
-       cdate=str(cdate)
-     else:
-       sys.exit('NDATE: Error - input cdate must be string or integer.  Exit!')
-   if not isinstance(hours, int):
-     if isinstance(hours, str):
-       hours=int(hours)
-     else:
-       sys.exit('NDATE: Error - input delta hour must be a string or integer.  Exit!')
-
-   indate=cdate.strip()
-   hh=indate[8:10]
-   yyyy=indate[0:4]
-   mm=indate[4:6]
-   dd=indate[6:8]
-   #set date/time field
-   parseme=(yyyy+' '+mm+' '+dd+' '+hh)
-   datetime_cdate=dateutil.parser.parse(parseme)
-   valid=datetime_cdate+dateutil.relativedelta.relativedelta(hours=+hours)
-   vyyyy=str(valid.year)
-   vm=str(valid.month).zfill(2)
-   vd=str(valid.day).zfill(2)
-   vh=str(valid.hour).zfill(2)
-   return vyyyy+vm+vd+vh
-
-
-def clear_plotables(ax,keep_ax_lst,fig):
-  #### - step to clear off old plottables but leave the map info - ####
-  if len(keep_ax_lst) == 0 :
-    print("clear_plotables WARNING keep_ax_lst has length 0. Clearing ALL plottables including map info!")
-  cur_ax_children = ax.get_children()[:]
-  if len(cur_ax_children) > 0:
-    for a in cur_ax_children:
-      if a not in keep_ax_lst:
-       # if the artist isn't part of the initial set up, remove it
-        a.remove()
-
-
-def compress_and_save(filename):
-  #### - compress and save the image - ####
-  ram = io.BytesIO()
-  plt.savefig(ram, format='png', bbox_inches='tight', dpi=150)
-  ram.seek(0)
-  im = Image.open(ram)
-  im2 = im.convert('RGB')#.convert('P', palette=Image.ADAPTIVE)
-  im2.save(filename, format='PNG')
-
-def cmap_t2m():
- # Create colormap for 2-m temperature
- # Modified version of the ncl_t2m colormap
-    r=np.array([255,128,0,  70, 51, 0,  255,0, 0,  51, 255,255,255,255,255,171,128,128,36,162,255])
-    g=np.array([0,  0,  0,  70, 102,162,255,92,128,185,255,214,153,102,0,  0,  0,  68, 36,162,255])
-    b=np.array([255,128,128,255,255,255,255,0, 0,  102,0,  112,0,  0,  0,  56, 0,  68, 36,162,255])
-    xsize=np.arange(np.size(r))
-    r = r/255.
-    g = g/255.
-    b = b/255.
-    red = []
-    green = []
-    blue = []
-    for i in range(len(xsize)):
-        xNorm=np.float(i)/(np.float(np.size(r))-1.0)
-        red.append([xNorm,r[i],r[i]])
-        green.append([xNorm,g[i],g[i]])
-        blue.append([xNorm,b[i],b[i]])
-    colorDict = {"red":red, "green":green, "blue":blue}
-    cmap_t2m_coltbl = matplotlib.colors.LinearSegmentedColormap('CMAP_T2M_COLTBL',colorDict)
-    return cmap_t2m_coltbl
-
-
-def cmap_q2m():
- # Create colormap for 2-m dew point temperature
-    r=np.array([255,179,96,128,0, 0,  51, 0,  0,  0,  133,51, 70, 0,  128,128,180])
-    g=np.array([255,179,96,128,92,128,153,155,155,255,162,102,70, 0,  0,  0,  0])
-    b=np.array([255,179,96,0,  0, 0,  102,155,255,255,255,255,255,128,255,128,128])
-    xsize=np.arange(np.size(r))
-    r = r/255.
-    g = g/255.
-    b = b/255.
-    red = []
-    green = []
-    blue = []
-    for i in range(len(xsize)):
-        xNorm=np.float(i)/(np.float(np.size(r))-1.0)
-        red.append([xNorm,r[i],r[i]])
-        green.append([xNorm,g[i],g[i]])
-        blue.append([xNorm,b[i],b[i]])
-    colorDict = {"red":red, "green":green, "blue":blue}
-    cmap_q2m_coltbl = matplotlib.colors.LinearSegmentedColormap('CMAP_Q2M_COLTBL',colorDict)
-    cmap_q2m_coltbl.set_over(color='deeppink')
-    return cmap_q2m_coltbl
-
-
-def rotate_wind(true_lat,lov_lon,earth_lons,uin,vin,proj,inverse=False):
-  #  Rotate winds from LCC relative to earth relative (or vice-versa if inverse==true)
-  #   This routine is vectorized and *should* work on any size 2D vg and ug arrays.
-  #   Program will quit if dimensions are too large.
-  #
-  # Input args:
-  #  true_lat = True latitidue for LCC projection (single value in degrees)
-  #  lov_lon  = The LOV value from grib (e.g. - -95.0) (single value in degrees)
-  #              Grib doc says: "Lov = orientation of the grid; i.e. the east longitude value of
-  #                              the meridian which is parallel to the Y-axis (or columns of the grid)
-  #                              along which latitude increases as the Y-coordinate increases (the
-  #                              orientation longitude may or may not appear on a particular grid).
-  #
-  #  earth_lons = Earth relative longitudes (can be an array, in degrees)
-  #  uin, vin     = Input winds to rotate
-  #
-  # Returns:
-  #  uout, vout = Output, rotated winds
-  #-----------------------------------------------------------------------------------------------------
-
-  # Get size and length of input u winds, if not 2d, raise an error
-  q=np.shape(uin)
-  ndims=len(q)
-  if ndims > 2:
-    # Raise error and quit!
-    raise SystemExit("Input winds for rotation have greater than 2 dimensions!")
-  if lov_lon > 0.: lov_lon=lov_lon-360.
-  dtr=np.pi/180.0             # Degrees to radians
-
-  if not isinstance(inverse, bool):
-    raise TypeError("**kwarg inverse must be of type bool.")
-
-  # Compute rotation constant which is also
-  # known as the Lambert cone constant.  In the case
-  # of a polar stereographic projection, this is one.
-  # See the following pdf for excellent documentation
-  # http://www.dtcenter.org/met/users/docs/write_ups/velocity.pdf
-  if proj.lower()=='lcc':
-    rotcon_p=np.sin(true_lat*dtr)
-  elif proj.lower() in ['stere','spstere', 'npstere']:
-    rotcon_p=1.0
-  else:
-    raise SystemExit("Unsupported map projection: "+proj.lower()+" for wind rotation.")
-
-  angles = rotcon_p*(earth_lons-lov_lon)*dtr
-  sinx2 = np.sin(angles)
-  cosx2 = np.cos(angles)
-
-  # Steps below are elementwise products, not matrix mutliplies
-  if inverse==False:
-    # Return the earth relative winds
-    uout = cosx2*uin+sinx2*vin
-    vout =-sinx2*uin+cosx2*vin
-  elif inverse==True:
-    # Return the grid relative winds
-    uout = cosx2*uin-sinx2*vin
-    vout = sinx2*uin+cosx2*vin
-
-  return uout,vout
-
-
-#-------------Start of script -------------------------#
-
-# Define required positional arguments
-parser = argparse.ArgumentParser()
-parser.add_argument("Cycle date/time in YYYYMMDDHH format")
-parser.add_argument("Starting forecast hour")
-parser.add_argument("Ending forecast hour")
-parser.add_argument("Forecast hour increment")
-parser.add_argument("Path to experiment directory")
-parser.add_argument("Path to base directory of cartopy shapefiles")
-parser.add_argument("Name of native domain used in forecast (and in constructing post file names)")
-args = parser.parse_args()
-
-# Read date/time, forecast hour, and directory paths from command line
-ymdh = str(sys.argv[1])
-ymd = ymdh[0:8]
-year = int(ymdh[0:4])
-month = int(ymdh[4:6])
-day = int(ymdh[6:8])
-hour = int(ymdh[8:10])
-cyc = str(hour).zfill(2)
-print(year, month, day, hour)
-
-# Define the range of forecast hours to create plots for
-start_fhr = int(sys.argv[2])
-end_fhr = int(sys.argv[3])
-increment_fhr = int(sys.argv[4])
-if (start_fhr == end_fhr) or (increment_fhr == 0):
-  fhours = [start_fhr]
-else:
-  num = int(((end_fhr - start_fhr) / increment_fhr) + 1)
-  fhours = np.linspace(start_fhr,end_fhr,num,dtype='int')
-print(fhours)
-
-EXPT_DIR = str(sys.argv[5])
-CARTOPY_DIR = str(sys.argv[6])
-POST_OUTPUT_DOMAIN_NAME = str(sys.argv[7]).lower()
-
-# Loop over forecast hours
-for fhr in fhours:
-  fhour = str(fhr).zfill(3)
-  print('Working on forecast hour '+fhour)
-  itime = ymdh
-  vtime = ndate(itime,int(fhr))
-
-# Define the location of the input file
-  data1 = pygrib.open(EXPT_DIR+'/'+ymdh+'/postprd/rrfs.t'+cyc+'z.prslev.f'+fhour+'.'+POST_OUTPUT_DOMAIN_NAME+'.grib2')
-
-# Get the lats and lons
-  grids = [data1]
-  lats = []
-  lons = []
-  lats_shift = []
-  lons_shift = []
-
-  for data in grids:
-    # Unshifted grid for contours and wind barbs
-    lat, lon = data[1].latlons()
-    lats.append(lat)
-    lons.append(lon)
-
-    # Shift grid for pcolormesh
-    lat1 = data[1]['latitudeOfFirstGridPointInDegrees']
-    lon1 = data[1]['longitudeOfFirstGridPointInDegrees']
-    try:
-        nx = data[1]['Nx']
-        ny = data[1]['Ny']
-    except:
-        nx = data[1]['Ni']
-        ny = data[1]['Nj']
-    dx = data[1]['DxInMetres']
-    dy = data[1]['DyInMetres']
-    pj = pyproj.Proj(data[1].projparams)
-    llcrnrx, llcrnry = pj(lon1,lat1)
-    llcrnrx = llcrnrx - (dx/2.)
-    llcrnry = llcrnry - (dy/2.)
-    x = llcrnrx + dx*np.arange(nx)
-    y = llcrnry + dy*np.arange(ny)
-    x,y = np.meshgrid(x,y)
-    lon, lat = pj(x, y, inverse=True)
-    lats_shift.append(lat)
-    lons_shift.append(lon)
-
-# Unshifted lat/lon arrays grabbed directly using latlons() method
-  lat = lats[0]
-  lon = lons[0]
-
-# Shifted lat/lon arrays for pcolormesh
-  lat_shift = lats_shift[0]
-  lon_shift = lons_shift[0]
-
-  Lat0 = data1[1]['LaDInDegrees']
-  Lon0 = data1[1]['LoVInDegrees']
-  print(Lat0)
-  print(Lon0)
-
-# Specify plotting domains
-# User can add domains here, just need to specify lat/lon information below
-# (if dom == 'conus' block)
-  domains=['conus']    # Other option is 'regional'
-
-###################################################
-# Read in all variables and calculate differences #
-###################################################
-  t1a = time.perf_counter()
-
-# Sea level pressure
-  slp = data1.select(name='Pressure reduced to MSL')[0].values * 0.01
-  slpsmooth = ndimage.filters.gaussian_filter(slp, 13.78)
-
-# 2-m temperature
-  tmp2m = data1.select(name='2 metre temperature')[0].values
-  tmp2m = (tmp2m - 273.15)*1.8 + 32.0
-
-# 2-m dew point temperature
-  dew2m = data1.select(name='2 metre dewpoint temperature')[0].values
-  dew2m = (dew2m - 273.15)*1.8 + 32.0
-
-# 10-m wind speed
-  uwind = data1.select(name='10 metre U wind component')[0].values * 1.94384
-  vwind = data1.select(name='10 metre V wind component')[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  uwind, vwind = rotate_wind(Lat0,Lon0,lon,uwind,vwind,'lcc',inverse=False)
-  wspd10m = np.sqrt(uwind**2 + vwind**2)
-
-# Surface-based CAPE
-  cape = data1.select(name='Convective available potential energy',typeOfLevel='surface')[0].values
-
-# Surface-based CIN
-  cin = data1.select(name='Convective inhibition',typeOfLevel='surface')[0].values
-
-# 500 mb height, wind, vorticity
-  z500 = data1.select(name='Geopotential height',level=500)[0].values * 0.1
-  z500 = ndimage.filters.gaussian_filter(z500, 6.89)
-  vort500 = data1.select(name='Absolute vorticity',level=500)[0].values * 100000
-  vort500 = ndimage.filters.gaussian_filter(vort500,1.7225)
-  vort500[vort500 > 1000] = 0	# Mask out undefined values on domain edge
-  u500 = data1.select(name='U component of wind',level=500)[0].values * 1.94384
-  v500 = data1.select(name='V component of wind',level=500)[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  u500, v500 = rotate_wind(Lat0,Lon0,lon,u500,v500,'lcc',inverse=False)
-
-# 250 mb winds
-  u250 = data1.select(name='U component of wind',level=250)[0].values * 1.94384
-  v250 = data1.select(name='V component of wind',level=250)[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  u250, v250 = rotate_wind(Lat0,Lon0,lon,u250,v250,'lcc',inverse=False)
-  wspd250 = np.sqrt(u250**2 + v250**2)
-
-# Total precipitation
-  qpf = data1.select(name='Total Precipitation',lengthOfTimeRange=fhr)[0].values * 0.0393701
-
-# Composite reflectivity
-  refc = data1.select(name='Maximum/Composite radar reflectivity')[0].values
-
-  if (fhr > 0):
-# Max/Min Hourly 2-5 km Updraft Helicity
-    maxuh25 = data1.select(stepType='max',parameterName="199",topLevel=5000,bottomLevel=2000)[0].values
-    minuh25 = data1.select(stepType='min',parameterName="200",topLevel=5000,bottomLevel=2000)[0].values
-    maxuh25[maxuh25 < 10] = 0
-    minuh25[minuh25 > -10] = 0
-    uh25 = maxuh25 + minuh25
-
-
-  t2a = time.perf_counter()
-  t3a = round(t2a-t1a, 3)
-  print(("%.3f seconds to read all messages") % t3a)
-
-
-########################################
-#    START PLOTTING FOR EACH DOMAIN    #
-########################################
-
-  def main():
-
-  # Number of processes must coincide with the number of domains to plot
-    #pool = multiprocessing.Pool(len(domains))
-    #pool.map(plot_all,domains)
-
-    # To avoid import multiprocessing recursively on MacOS etc.
-    # Anyway since we only have one domain for SRW application
-    for dom in domains:
-        plot_all(dom)
-
-  def plot_all(dom):
-
-    t1dom = time.perf_counter()
-
-  # Map corners for each domain
-    if dom == 'conus':
-      llcrnrlon = -120.5
-      llcrnrlat = 21.0
-      urcrnrlon = -64.5
-      urcrnrlat = 49.0
-      lat_0 = 35.4
-      lon_0 = -97.6
-      extent=[llcrnrlon-3,urcrnrlon-6,llcrnrlat-1,urcrnrlat+2]
-    elif dom == 'regional':
-      llcrnrlon = np.min(lon)
-      llcrnrlat = np.min(lat)
-      urcrnrlon = np.max(lon)
-      urcrnrlat = np.max(lat)
-      lat_0 = Lat0
-      lon_0 = Lon0
-      extent=[llcrnrlon,urcrnrlon,llcrnrlat-1,urcrnrlat]
-
-
-  # create figure and axes instances
-    fig = plt.figure(figsize=(10,10))
-    ax1 = fig.add_axes([0.1,0.1,0.8,0.8])
-
-  # Define where Cartopy Maps are located
-    cartopy.config['data_dir'] = CARTOPY_DIR
-
-    back_res='50m'
-    back_img='on'
-
-  # set up the map background with cartopy
-    myproj=ccrs.LambertConformal(central_longitude=lon_0, central_latitude=lat_0, false_easting=0.0,
-                            false_northing=0.0, secant_latitudes=None, standard_parallels=None,
-                            globe=None)
-    ax = plt.axes(projection=myproj)
-    ax.set_extent(extent)
-
-    fline_wd = 0.5  # line width
-    falpha = 0.3    # transparency
-
-  # natural_earth
-#  land=cfeature.NaturalEarthFeature('physical','land',back_res,
-#                    edgecolor='face',facecolor=cfeature.COLORS['land'],
-#                    alpha=falpha)
-    lakes=cfeature.NaturalEarthFeature('physical','lakes',back_res,
-                      edgecolor='blue',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-    coastline=cfeature.NaturalEarthFeature('physical','coastline',
-                      back_res,edgecolor='blue',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-    states=cfeature.NaturalEarthFeature('cultural','admin_1_states_provinces',
-                      back_res,edgecolor='black',facecolor='none',
-                      linewidth=fline_wd,linestyle=':',alpha=falpha)
-    borders=cfeature.NaturalEarthFeature('cultural','admin_0_countries',
-                      back_res,edgecolor='red',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-
-  # All lat lons are earth relative, so setup the associated projection correct for that data
-    transform = ccrs.PlateCarree()
-
-  # high-resolution background images
-    if back_img=='on':
-       img = plt.imread(CARTOPY_DIR+'/raster_files/NE1_50M_SR_W.tif')
-       ax.imshow(img, origin='upper', transform=transform)
-
-#  ax.add_feature(land)
-    ax.add_feature(lakes)
-    ax.add_feature(states)
-    ax.add_feature(borders)
-    ax.add_feature(coastline)
-
-  # Map/figure has been set up here, save axes instances for use again later
-    keep_ax_lst = ax.get_children()[:]
-
-
-################################
-  # Plot SLP
-################################
-    t1 = time.perf_counter()
-    print(('Working on slp for '+dom))
-
-    units = 'mb'
-    clevs = [976,980,984,988,992,996,1000,1004,1008,1012,1016,1020,1024,1028,1032,1036,1040,1044,1048,1052]
-    clevsdif = [-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12]
-    cm = plt.cm.Spectral_r
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs1_a = plt.pcolormesh(lon_shift,lat_shift,slp,transform=transform,cmap=cm,norm=norm)
-    cbar1 = plt.colorbar(cs1_a,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    cs1_b = plt.contour(lon_shift,lat_shift,slpsmooth,np.arange(940,1060,4),colors='black',linewidths=1.25,transform=transform)
-    plt.clabel(cs1_b,np.arange(940,1060,4),inline=1,fmt='%d',fontsize=8)
-    ax.text(.5,1.03,'FV3-LAM SLP ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/slp_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot slp for: '+dom) % t3)
-
-
-#################################
-  # Plot 2-m T
-#################################
-    t1 = time.perf_counter()
-    print(('Working on t2m for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = '\xb0''F'
-    clevs = np.linspace(-16,134,51)
-    cm = plt.cm.Spectral_r #cmap_t2m()
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,tmp2m,transform=transform,cmap=cm,norm=norm)
-    cs_1.cmap.set_under('white')
-    cs_1.cmap.set_over('white')
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=[-16,-4,8,20,32,44,56,68,80,92,104,116,128],extend='both')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    ax.text(.5,1.03,'FV3-LAM 2-m Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/2mt_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 2mt for: '+dom) % t3)
-
-
-#################################
-  # Plot 2-m Dew Point
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 2mdew for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = '\xb0''F'
-    clevs = np.linspace(-5,80,35)
-    cm = cmap_q2m()
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,dew2m,transform=transform,cmap=cm,norm=norm)
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    ax.text(.5,1.03,'FV3-LAM 2-m Dew Point Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/2mdew_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 2mdew for: '+dom) % t3)
-
-
-#################################
-  # Plot 10-m WSPD
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 10mwspd for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = 'kts'
-  # Places a wind barb every ~180 km, optimized for CONUS domain
-    skip = round(177.28*(dx/1000.)**-.97)
-    print('skipping every '+str(skip)+' grid points to plot')
-    barblength = 4
-
-    clevs = [5,10,15,20,25,30,35,40,45,50,55,60]
-    colorlist = ['turquoise','dodgerblue','blue','#FFF68F','#E3CF57','peru','brown','crimson','red','fuchsia','DarkViolet']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,wspd10m,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    plt.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],uwind[::skip,::skip],vwind[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax.text(.5,1.03,'FV3-LAM 10-m Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/10mwind_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 10mwspd for: '+dom) % t3)
-
-
-#################################
-  # Plot Surface-Based CAPE/CIN
-#################################
-    t1 = time.perf_counter()
-    print(('Working on surface-based CAPE/CIN for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = 'J/kg'
-    clevs = [100,250,500,1000,1500,2000,2500,3000,3500,4000,4500,5000]
-    clevs2 = [-2000,-500,-250,-100,-25]
-    colorlist = ['blue','dodgerblue','cyan','mediumspringgreen','#FAFAD2','#EEEE00','#EEC900','darkorange','crimson','darkred','darkviolet']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,cape,transform=transform,cmap=cm,vmin=100,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    cs_1b = plt.contourf(lon_shift,lat_shift,cin,clevs2,colors='none',hatches=['**','++','////','..'],transform=transform)
-    ax.text(.5,1.05,'FV3-LAM Surface-Based CAPE (shaded) and CIN (hatched) ('+units+') \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/sfcape_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot surface-based CAPE/CIN for: '+dom) % t3)
-
-
-#################################
-  # Plot 500 mb HGT/WIND/VORT
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 500 mb Hgt/Wind/Vort for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = 'x10${^5}$ s${^{-1}}$'
-    skip = round(177.28*(dx/1000.)**-.97)
-    barblength = 4
-
-    vortlevs = [16,20,24,28,32,36,40]
-    colorlist = ['yellow','gold','goldenrod','orange','orangered','red']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(vortlevs, cm.N)
-
-    cs1_a = plt.pcolormesh(lon_shift,lat_shift,vort500,transform=transform,cmap=cm,norm=norm)
-    cs1_a.cmap.set_under('white')
-    cs1_a.cmap.set_over('darkred')
-    cbar1 = plt.colorbar(cs1_a,orientation='horizontal',pad=0.05,shrink=0.6,ticks=vortlevs,extend='both')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    plt.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],u500[::skip,::skip],v500[::skip,::skip],length=barblength,linewidth=0.5,color='steelblue',transform=transform)
-    cs1_b = plt.contour(lon_shift,lat_shift,z500,np.arange(486,600,6),colors='black',linewidths=1,transform=transform)
-    plt.clabel(cs1_b,np.arange(486,600,6),inline_spacing=1,fmt='%d',fontsize=8)
-    ax.text(.5,1.03,'FV3-LAM 500 mb Heights (dam), Winds (kts), and $\zeta$ ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/500_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 500 mb Hgt/Wind/Vort for: '+dom) % t3)
-
-
-#################################
-  # Plot 250 mb WIND
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 250 mb WIND for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = 'kts'
-    skip = round(177.28*(dx/1000.)**-.97)
-
-    barblength = 4
-
-    clevs = [50,60,70,80,90,100,110,120,130,140,150]
-    colorlist = ['turquoise','deepskyblue','dodgerblue','#1874CD','blue','beige','khaki','peru','brown','crimson']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,wspd250,transform=transform,cmap=cm,vmin=50,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('red')
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    plt.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],u250[::skip,::skip],v250[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax.text(.5,1.03,'FV3-LAM 250 mb Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/250wind_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 250 mb WIND for: '+dom) % t3)
-
-#################################
-  # Plot Total QPF
-#################################
-    if (fhr > 0):		# Do not make total QPF plot for forecast hour 0
-      t1 = time.perf_counter()
-      print(('Working on total qpf for '+dom))
-
-    # Clear off old plottables but keep all the map info
-      cbar1.remove()
-      clear_plotables(ax,keep_ax_lst,fig)
-
-      units = 'in'
-      clevs = [0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,7,10,15,20]
-      clevsdif = [-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3]
-      colorlist = ['chartreuse','limegreen','green','blue','dodgerblue','deepskyblue','cyan','mediumpurple','mediumorchid','darkmagenta','darkred','crimson','orangered','darkorange','goldenrod','gold','yellow']
-      cm = matplotlib.colors.ListedColormap(colorlist)
-      norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-      cs_1 = plt.pcolormesh(lon_shift,lat_shift,qpf,transform=transform,cmap=cm,vmin=0.01,norm=norm)
-      cs_1.cmap.set_under('white',alpha=0.)
-      cs_1.cmap.set_over('pink')
-      cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-      cbar1.set_label(units,fontsize=8)
-      cbar1.ax.set_xticklabels(clevs)
-      cbar1.ax.tick_params(labelsize=8)
-      ax.text(.5,1.03,'FV3-LAM '+fhour+'-hr Accumulated Precipitation ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/qpf_'+dom+'_f'+fhour+'.png')
-      t2 = time.perf_counter()
-      t3 = round(t2-t1, 3)
-      print(('%.3f seconds to plot total qpf for: '+dom) % t3)
-
-#################################
-  # Plot composite reflectivity
-#################################
-    t1 = time.perf_counter()
-    print(('Working on composite reflectivity for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    clear_plotables(ax,keep_ax_lst,fig)
-
-    units = 'dBZ'
-    clevs = np.linspace(5,70,14)
-    clevsdif = [20,1000]
-    colorlist = ['turquoise','dodgerblue','mediumblue','lime','limegreen','green','#EEEE00','#EEC900','darkorange','red','firebrick','darkred','fuchsia']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = plt.pcolormesh(lon_shift,lat_shift,refc,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=8)
-    cbar1.ax.tick_params(labelsize=8)
-    ax.text(.5,1.03,'FV3-LAM Composite Reflectivity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/refc_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot composite reflectivity for: '+dom) % t3)
-
-
-#################################
-  # Plot Max/Min Hourly 2-5 km UH
-#################################
-    if (fhr > 0):		# Do not make max/min hourly 2-5 km UH plot for forecast hour 0
-      t1 = time.perf_counter()
-      print(('Working on Max/Min Hourly 2-5 km UH for '+dom))
-
-    # Clear off old plottables but keep all the map info
-      cbar1.remove()
-      clear_plotables(ax,keep_ax_lst,fig)
-
-      units = 'm${^2}$ s$^{-2}$'
-      clevs = [-150,-100,-75,-50,-25,-10,0,10,25,50,75,100,150,200,250,300]
-#   alternative colormap for just max UH if you don't want to plot the min UH too
-#     colorlist = ['white','skyblue','mediumblue','green','orchid','firebrick','#EEC900','DarkViolet']
-      colorlist = ['blue','#1874CD','dodgerblue','deepskyblue','turquoise','#E5E5E5','#E5E5E5','#EEEE00','#EEC900','darkorange','orangered','red','firebrick','mediumvioletred','darkviolet']
-      cm = matplotlib.colors.ListedColormap(colorlist)
-      norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-      cs_1 = plt.pcolormesh(lon_shift,lat_shift,uh25,transform=transform,cmap=cm,norm=norm)
-      cs_1.cmap.set_under('darkblue')
-      cs_1.cmap.set_over('black')
-      cbar1 = plt.colorbar(cs_1,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-      cbar1.set_label(units,fontsize=8)
-      cbar1.ax.tick_params(labelsize=8)
-      ax.text(.5,1.03,'FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      compress_and_save(EXPT_DIR+'/'+ymdh+'/postprd/uh25_'+dom+'_f'+fhour+'.png')
-      t2 = time.perf_counter()
-      t3 = round(t2-t1, 3)
-      print(('%.3f seconds to plot Max/Min Hourly 2-5 km UH for: '+dom) % t3)
-
-
-######################################################
-
-    t3dom = round(t2-t1dom, 3)
-    print(("%.3f seconds to plot all variables for forecast hour "+fhour) % t3dom)
-    plt.clf()
-
-######################################################
-
-  main()
-
diff --git a/ush/Python/plot_allvars_diff.py b/ush/Python/plot_allvars_diff.py
deleted file mode 100755
index 12fcdbd949..0000000000
--- a/ush/Python/plot_allvars_diff.py
+++ /dev/null
@@ -1,1101 +0,0 @@
-################################################################################
-####  Python Script Documentation Block
-#
-# Script name:       	plot_allvars_diff.py
-# Script description:  	Generates difference plots from FV3-LAM post processed
-#                       grib2 output over the CONUS
-#
-# Authors:  Ben Blake		Org: NOAA/NWS/NCEP/EMC		Date: 2020-08-24
-#           David Wright 	Org: University of Michigan
-#
-# Instructions:		Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour
-#                       3. Ending forecast hour
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR_1: Experiment 1 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_1/YYYYMMDDHH/postprd/
-#                       6. EXPT_DIR_2: Experiment 2 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_2/YYYYMMDDHH/postprd/
-#                       7. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -Shapefiles cannot be directly downloaded to NOAA
-#                            machines from the internet, so shapefiles need to
-#                            be downloaded if geopolitical boundaries are
-#                            desired on the maps.
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                          -More information regarding files needed to setup
-#                            display maps in Cartopy, see SRW App Users' Guide
-#                       8. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecasts and in constructing the names 
-#                          of the post output files.  This must be the same 
-#                          for both forecasts.
-#
-#           		To create difference plots for two forecasts on the
-#                       RRFS_CONUS_25km grid for hours 20-24 from the 5/7
-#                       00Z cycle with hourly output:
-#                          python plot_allvars_diff.py 2020050700  20 24 1 \
-#                          /path/to/expt_dir_1 /path/to/expt_dir_2 \
-#                          /path/to/base/cartopy/maps RRFS_CONUS_25km
-#
-#                       The variable domains in this script can be set to either
-#                         'conus' for a CONUS map or 'regional' where the map
-#                         is defined from variables in the grib2 files
-#
-################################################################################
-
-#-------------Import modules --------------------------#
-import pygrib
-import cartopy.crs as ccrs
-from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
-import cartopy.feature as cfeature
-import matplotlib
-matplotlib.use('Agg')
-import io
-import matplotlib.pyplot as plt
-import dateutil.relativedelta, dateutil.parser
-from PIL import Image
-from matplotlib.gridspec import GridSpec
-import numpy as np
-import time,os,sys,multiprocessing
-import multiprocessing.pool
-from scipy import ndimage
-import pyproj
-import argparse
-import cartopy
-
-#--------------Define some functions ------------------#
-
-def ndate(cdate,hours):
-   if not isinstance(cdate, str):
-     if isinstance(cdate, int):
-       cdate=str(cdate)
-     else:
-       sys.exit('NDATE: Error - input cdate must be string or integer.  Exit!')
-   if not isinstance(hours, int):
-     if isinstance(hours, str):
-       hours=int(hours)
-     else:
-       sys.exit('NDATE: Error - input delta hour must be a string or integer.  Exit!')
-
-   indate=cdate.strip()
-   hh=indate[8:10]
-   yyyy=indate[0:4]
-   mm=indate[4:6]
-   dd=indate[6:8]
-   #set date/time field
-   parseme=(yyyy+' '+mm+' '+dd+' '+hh)
-   datetime_cdate=dateutil.parser.parse(parseme)
-   valid=datetime_cdate+dateutil.relativedelta.relativedelta(hours=+hours)
-   vyyyy=str(valid.year)
-   vm=str(valid.month).zfill(2)
-   vd=str(valid.day).zfill(2)
-   vh=str(valid.hour).zfill(2)
-   return vyyyy+vm+vd+vh
-
-
-def clear_plotables(ax,keep_ax_lst,fig):
-  #### - step to clear off old plottables but leave the map info - ####
-  if len(keep_ax_lst) == 0 :
-    print("clear_plotables WARNING keep_ax_lst has length 0. Clearing ALL plottables including map info!")
-  cur_ax_children = ax.get_children()[:]
-  if len(cur_ax_children) > 0:
-    for a in cur_ax_children:
-      if a not in keep_ax_lst:
-       # if the artist isn't part of the initial set up, remove it
-        a.remove()
-
-
-def compress_and_save(filename):
-  #### - compress and save the image - ####
-  ram = io.BytesIO()
-  plt.savefig(ram, format='png', bbox_inches='tight', dpi=150)
-  ram.seek(0)
-  im = Image.open(ram)
-  im2 = im.convert('RGB')#.convert('P', palette=Image.ADAPTIVE)
-  im2.save(filename, format='PNG')
-
-def cmap_t2m():
- # Create colormap for 2-m temperature
- # Modified version of the ncl_t2m colormap
-    r=np.array([255,128,0,  70, 51, 0,  255,0, 0,  51, 255,255,255,255,255,171,128,128,36,162,255])
-    g=np.array([0,  0,  0,  70, 102,162,255,92,128,185,255,214,153,102,0,  0,  0,  68, 36,162,255])
-    b=np.array([255,128,128,255,255,255,255,0, 0,  102,0,  112,0,  0,  0,  56, 0,  68, 36,162,255])
-    xsize=np.arange(np.size(r))
-    r = r/255.
-    g = g/255.
-    b = b/255.
-    red = []
-    green = []
-    blue = []
-    for i in range(len(xsize)):
-        xNorm=np.float(i)/(np.float(np.size(r))-1.0)
-        red.append([xNorm,r[i],r[i]])
-        green.append([xNorm,g[i],g[i]])
-        blue.append([xNorm,b[i],b[i]])
-    colorDict = {"red":red, "green":green, "blue":blue}
-    cmap_t2m_coltbl = matplotlib.colors.LinearSegmentedColormap('CMAP_T2M_COLTBL',colorDict)
-    return cmap_t2m_coltbl
-
-
-def cmap_q2m():
- # Create colormap for 2-m dew point temperature
-    r=np.array([255,179,96,128,0, 0,  51, 0,  0,  0,  133,51, 70, 0,  128,128,180])
-    g=np.array([255,179,96,128,92,128,153,155,155,255,162,102,70, 0,  0,  0,  0])
-    b=np.array([255,179,96,0,  0, 0,  102,155,255,255,255,255,255,128,255,128,128])
-    xsize=np.arange(np.size(r))
-    r = r/255.
-    g = g/255.
-    b = b/255.
-    red = []
-    green = []
-    blue = []
-    for i in range(len(xsize)):
-        xNorm=np.float(i)/(np.float(np.size(r))-1.0)
-        red.append([xNorm,r[i],r[i]])
-        green.append([xNorm,g[i],g[i]])
-        blue.append([xNorm,b[i],b[i]])
-    colorDict = {"red":red, "green":green, "blue":blue}
-    cmap_q2m_coltbl = matplotlib.colors.LinearSegmentedColormap('CMAP_Q2M_COLTBL',colorDict)
-    cmap_q2m_coltbl.set_over(color='deeppink')
-    return cmap_q2m_coltbl
-
-
-def rotate_wind(true_lat,lov_lon,earth_lons,uin,vin,proj,inverse=False):
-  #  Rotate winds from LCC relative to earth relative (or vice-versa if inverse==true)
-  #   This routine is vectorized and *should* work on any size 2D vg and ug arrays.
-  #   Program will quit if dimensions are too large.
-  #
-  # Input args:
-  #  true_lat = True latitidue for LCC projection (single value in degrees)
-  #  lov_lon  = The LOV value from grib (e.g. - -95.0) (single value in degrees)
-  #              Grib doc says: "Lov = orientation of the grid; i.e. the east longitude value of
-  #                              the meridian which is parallel to the Y-axis (or columns of the grid)
-  #                              along which latitude increases as the Y-coordinate increases (the
-  #                              orientation longitude may or may not appear on a particular grid).
-  #
-  #  earth_lons = Earth relative longitudes (can be an array, in degrees)
-  #  uin, vin     = Input winds to rotate
-  #
-  # Returns:
-  #  uout, vout = Output, rotated winds
-  #-----------------------------------------------------------------------------------------------------
-
-  # Get size and length of input u winds, if not 2d, raise an error
-  q=np.shape(uin)
-  ndims=len(q)
-  if ndims > 2:
-    # Raise error and quit!
-    raise SystemExit("Input winds for rotation have greater than 2 dimensions!")
-  if lov_lon > 0.: lov_lon=lov_lon-360.
-  dtr=np.pi/180.0             # Degrees to radians
-
-  if not isinstance(inverse, bool):
-    raise TypeError("**kwarg inverse must be of type bool.")
-
-  # Compute rotation constant which is also
-  # known as the Lambert cone constant.  In the case
-  # of a polar stereographic projection, this is one.
-  # See the following pdf for excellent documentation
-  # http://www.dtcenter.org/met/users/docs/write_ups/velocity.pdf
-  if proj.lower()=='lcc':
-    rotcon_p=np.sin(true_lat*dtr)
-  elif proj.lower() in ['stere','spstere', 'npstere']:
-    rotcon_p=1.0
-  else:
-    raise SystemExit("Unsupported map projection: "+proj.lower()+" for wind rotation.")
-
-  angles = rotcon_p*(earth_lons-lov_lon)*dtr
-  sinx2 = np.sin(angles)
-  cosx2 = np.cos(angles)
-
-  # Steps below are elementwise products, not matrix mutliplies
-  if inverse==False:
-    # Return the earth relative winds
-    uout = cosx2*uin+sinx2*vin
-    vout =-sinx2*uin+cosx2*vin
-  elif inverse==True:
-    # Return the grid relative winds
-    uout = cosx2*uin-sinx2*vin
-    vout = sinx2*uin+cosx2*vin
-
-  return uout,vout
-
-
-#-------------Start of script -------------------------#
-
-# Define required positional arguments
-parser = argparse.ArgumentParser()
-parser.add_argument("Cycle date/time in YYYYMMDDHH format")
-parser.add_argument("Starting forecast hour")
-parser.add_argument("Ending forecast hour")
-parser.add_argument("Forecast hour increment")
-parser.add_argument("Path to experiment 1 directory")
-parser.add_argument("Path to experiment 2 directory")
-parser.add_argument("Path to base directory of cartopy shapefiles")
-parser.add_argument("Name of native domain used in forecasts (and in constructing post file names)")
-args = parser.parse_args()
-
-# Read date/time, forecast hour, and directory paths from command line
-ymdh = str(sys.argv[1])
-ymd = ymdh[0:8]
-year = int(ymdh[0:4])
-month = int(ymdh[4:6])
-day = int(ymdh[6:8])
-hour = int(ymdh[8:10])
-cyc = str(hour).zfill(2)
-print(year, month, day, hour)
-
-# Define the range of forecast hours to create plots for
-start_fhr = int(sys.argv[2])
-end_fhr = int(sys.argv[3])
-increment_fhr = int(sys.argv[4])
-if (start_fhr == end_fhr) or (increment_fhr == 0):
-  fhours = [start_fhr]
-else:
-  num = int(((end_fhr - start_fhr) / increment_fhr) + 1)
-  fhours = np.linspace(start_fhr,end_fhr,num,dtype='int')
-print(fhours)
-
-EXPT_DIR_1 = str(sys.argv[5])
-EXPT_DIR_2 = str(sys.argv[6])
-CARTOPY_DIR = str(sys.argv[7])
-POST_OUTPUT_DOMAIN_NAME = str(sys.argv[8]).lower()
-
-# Loop over forecast hours
-for fhr in fhours:
-  fhour = str(fhr).zfill(3)
-  print('Working on forecast hour '+fhour)
-  itime = ymdh
-  vtime = ndate(itime,int(fhr))
-
-
-# Define the location of the input files
-  data1 = pygrib.open(EXPT_DIR_1+'/'+ymdh+'/postprd/rrfs.t'+cyc+'z.prslev.f'+fhour+'.'+POST_OUTPUT_DOMAIN_NAME+'.grib2')
-  data2 = pygrib.open(EXPT_DIR_2+'/'+ymdh+'/postprd/rrfs.t'+cyc+'z.prslev.f'+fhour+'.'+POST_OUTPUT_DOMAIN_NAME+'.grib2')
-
-# Get the lats and lons
-  grids = [data1, data2]
-  lats = []
-  lons = []
-  lats_shift = []
-  lons_shift = []
-
-  for data in grids:
-    # Unshifted grid for contours and wind barbs
-    lat, lon = data[1].latlons()
-    lats.append(lat)
-    lons.append(lon)
-
-    # Shift grid for pcolormesh
-    lat1 = data[1]['latitudeOfFirstGridPointInDegrees']
-    lon1 = data[1]['longitudeOfFirstGridPointInDegrees']
-    try:
-        nx = data[1]['Nx']
-        ny = data[1]['Ny']
-    except:
-        nx = data[1]['Ni']
-        ny = data[1]['Nj']
-    dx = data[1]['DxInMetres']
-    dy = data[1]['DyInMetres']
-    pj = pyproj.Proj(data[1].projparams)
-    llcrnrx, llcrnry = pj(lon1,lat1)
-    llcrnrx = llcrnrx - (dx/2.)
-    llcrnry = llcrnry - (dy/2.)
-    x = llcrnrx + dx*np.arange(nx)
-    y = llcrnry + dy*np.arange(ny)
-    x,y = np.meshgrid(x,y)
-    lon, lat = pj(x, y, inverse=True)
-    lats_shift.append(lat)
-    lons_shift.append(lon)
-
-# Unshifted lat/lon arrays grabbed directly using latlons() method
-  lat = lats[0]
-  lon = lons[0]
-  lat2 = lats[1]
-  lon2 = lons[1]
-
-# Shifted lat/lon arrays for pcolormesh
-  lat_shift = lats_shift[0]
-  lon_shift = lons_shift[0]
-  lat2_shift = lats_shift[1]
-  lon2_shift = lons_shift[1]
-
-  Lat0 = data1[1]['LaDInDegrees']
-  Lon0 = data1[1]['LoVInDegrees']
-  print(Lat0)
-  print(Lon0)
-
-# Specify plotting domains
-# User can add domains here, just need to specify lat/lon information below
-# (if dom == 'conus' block)
-  domains=['conus']    # Other option is 'regional'
-
-###################################################
-# Read in all variables and calculate differences #
-###################################################
-  t1a = time.perf_counter()
-
-# Sea level pressure
-  slp_1 = data1.select(name='Pressure reduced to MSL')[0].values * 0.01
-  slpsmooth_1 = ndimage.filters.gaussian_filter(slp_1, 13.78)
-  slp_2 = data2.select(name='Pressure reduced to MSL')[0].values * 0.01
-  slpsmooth_2 = ndimage.filters.gaussian_filter(slp_2, 13.78)
-  slp_diff = slp_2 - slp_1
-
-# 2-m temperature
-  tmp2m_1 = data1.select(name='2 metre temperature')[0].values
-  tmp2m_1 = (tmp2m_1 - 273.15)*1.8 + 32.0
-  tmp2m_2 = data2.select(name='2 metre temperature')[0].values
-  tmp2m_2 = (tmp2m_2 - 273.15)*1.8 + 32.0
-  tmp2m_diff = tmp2m_2 - tmp2m_1
-
-# 2-m dew point temperature
-  dew2m_1 = data1.select(name='2 metre dewpoint temperature')[0].values
-  dew2m_1 = (dew2m_1 - 273.15)*1.8 + 32.0
-  dew2m_2 = data2.select(name='2 metre dewpoint temperature')[0].values
-  dew2m_2 = (dew2m_2 - 273.15)*1.8 + 32.0
-  dew2m_diff = dew2m_2 - dew2m_1
-
-# 10-m wind speed
-  uwind_1 = data1.select(name='10 metre U wind component')[0].values * 1.94384
-  vwind_1 = data1.select(name='10 metre V wind component')[0].values * 1.94384
-  uwind_2 = data2.select(name='10 metre U wind component')[0].values * 1.94384
-  vwind_2 = data2.select(name='10 metre V wind component')[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  uwind_1, vwind_1 = rotate_wind(Lat0,Lon0,lon,uwind_1,vwind_1,'lcc',inverse=False)
-  uwind_2, vwind_2 = rotate_wind(Lat0,Lon0,lon2,uwind_2,vwind_2,'lcc',inverse=False)
-  wspd10m_1 = np.sqrt(uwind_1**2 + vwind_1**2)
-  wspd10m_2 = np.sqrt(uwind_2**2 + vwind_2**2)
-  wspd10m_diff = wspd10m_2 - wspd10m_1
-
-# Surface-based CAPE
-  cape_1 = data1.select(name='Convective available potential energy',typeOfLevel='surface')[0].values
-  cape_2 = data2.select(name='Convective available potential energy',typeOfLevel='surface')[0].values
-  cape_diff = cape_2 - cape_1
-
-# Surface-based CIN
-  cin_1 = data1.select(name='Convective inhibition',typeOfLevel='surface')[0].values
-  cin_2 = data2.select(name='Convective inhibition',typeOfLevel='surface')[0].values
-  cin_diff = cin_2 - cin_1
-
-# 500 mb height, wind, vorticity
-  z500_1 = data1.select(name='Geopotential height',level=500)[0].values * 0.1
-  z500_1 = ndimage.filters.gaussian_filter(z500_1, 6.89)
-  z500_2 = data2.select(name='Geopotential height',level=500)[0].values * 0.1
-  z500_2 = ndimage.filters.gaussian_filter(z500_2, 6.89)
-  z500_diff = z500_2 - z500_1
-  vort500_1 = data1.select(name='Absolute vorticity',level=500)[0].values * 100000
-  vort500_1 = ndimage.filters.gaussian_filter(vort500_1,1.7225)
-  vort500_1[vort500_1 > 1000] = 0	# Mask out undefined values on domain edge
-  vort500_2 = data2.select(name='Absolute vorticity',level=500)[0].values * 100000
-  vort500_2 = ndimage.filters.gaussian_filter(vort500_2,1.7225)
-  vort500_2[vort500_2 > 1000] = 0	# Mask out undefined values on domain edge
-  u500_1 = data1.select(name='U component of wind',level=500)[0].values * 1.94384
-  u500_2 = data2.select(name='U component of wind',level=500)[0].values * 1.94384
-  v500_1 = data1.select(name='V component of wind',level=500)[0].values * 1.94384
-  v500_2 = data2.select(name='V component of wind',level=500)[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  u500_1, v500_1 = rotate_wind(Lat0,Lon0,lon,u500_1,v500_1,'lcc',inverse=False)
-  u500_2, v500_2 = rotate_wind(Lat0,Lon0,lon2,u500_2,v500_2,'lcc',inverse=False)
-
-# 250 mb winds
-  u250_1 = data1.select(name='U component of wind',level=250)[0].values * 1.94384
-  u250_2 = data2.select(name='U component of wind',level=250)[0].values * 1.94384
-  v250_1 = data1.select(name='V component of wind',level=250)[0].values * 1.94384
-  v250_2 = data2.select(name='V component of wind',level=250)[0].values * 1.94384
-# Rotate winds from grid relative to Earth relative
-  u250_1, v250_1 = rotate_wind(Lat0,Lon0,lon,u250_1,v250_1,'lcc',inverse=False)
-  u250_2, v250_2 = rotate_wind(Lat0,Lon0,lon2,u250_2,v250_2,'lcc',inverse=False)
-  wspd250_1 = np.sqrt(u250_1**2 + v250_1**2)
-  wspd250_2 = np.sqrt(u250_2**2 + v250_2**2)
-  wspd250_diff = wspd250_2 - wspd250_1
-
-# Total precipitation
-  qpf_1 = data1.select(name='Total Precipitation',lengthOfTimeRange=fhr)[0].values * 0.0393701
-  qpf_2 = data2.select(name='Total Precipitation',lengthOfTimeRange=fhr)[0].values * 0.0393701
-  qpf_diff = qpf_2 - qpf_1
-
-# Composite reflectivity
-  refc_1 = data1.select(name='Maximum/Composite radar reflectivity')[0].values
-  refc_2 = data2.select(name='Maximum/Composite radar reflectivity')[0].values
-
-  if (fhr > 0):
-# Max/Min Hourly 2-5 km Updraft Helicity
-    maxuh25_1 = data1.select(stepType='max',parameterName="199",topLevel=5000,bottomLevel=2000)[0].values
-    maxuh25_2 = data2.select(stepType='max',parameterName="199",topLevel=5000,bottomLevel=2000)[0].values
-    minuh25_1 = data1.select(stepType='min',parameterName="200",topLevel=5000,bottomLevel=2000)[0].values
-    minuh25_2 = data2.select(stepType='min',parameterName="200",topLevel=5000,bottomLevel=2000)[0].values
-    maxuh25_1[maxuh25_1 < 10] = 0
-    maxuh25_2[maxuh25_2 < 10] = 0
-    minuh25_1[minuh25_1 > -10] = 0
-    minuh25_2[minuh25_2 > -10] = 0
-    uh25_1 = maxuh25_1 + minuh25_1
-    uh25_2 = maxuh25_2 + minuh25_2
-    uh25_diff = uh25_2 - uh25_1
-
-
-  t2a = time.perf_counter()
-  t3a = round(t2a-t1a, 3)
-  print(("%.3f seconds to read all messages") % t3a)
-
-
-########################################
-#    START PLOTTING FOR EACH DOMAIN    #
-########################################
-
-  def main():
-
-    # Number of processes must coincide with the number of domains to plot
-    #pool = multiprocessing.Pool(len(domains))
-    #pool.map(plot_all,domains)
-
-    # To avoid import multiprocessing recursively on MacOS/Windows etc.
-    # Anyway since we only have one domain for SRW application
-    for dom in domains:
-        plot_all(dom)
-
-  def plot_all(dom):
-
-    t1dom = time.perf_counter()
-
-  # Map corners for each domain
-    if dom == 'conus':
-      llcrnrlon = -120.5
-      llcrnrlat = 21.0
-      urcrnrlon = -64.5
-      urcrnrlat = 49.0
-      lat_0 = 35.4
-      lon_0 = -97.6
-      extent=[llcrnrlon-3,urcrnrlon-6,llcrnrlat-1,urcrnrlat+2]
-    elif dom == 'regional':
-      llcrnrlon = np.min(lon)
-      llcrnrlat = np.min(lat)
-      urcrnrlon = np.max(lon)
-      urcrnrlat = np.max(lat)
-      lat_0 = Lat0
-      lon_0 = Lon0
-      extent=[llcrnrlon,urcrnrlon,llcrnrlat-1,urcrnrlat]
-
-
-  # create figure and axes instances
-    fig = plt.figure(figsize=(10,10))
-    gs = GridSpec(9,9,wspace=0.0,hspace=0.0)
-
-  # Define where Cartopy Maps are located
-    cartopy.config['data_dir'] = CARTOPY_DIR
-
-    back_res='50m'
-    back_img='on'
-
-  # set up the map background with cartopy
-    myproj=ccrs.LambertConformal(central_longitude=lon_0, central_latitude=lat_0, false_easting=0.0,
-                            false_northing=0.0, secant_latitudes=None, standard_parallels=None,
-                            globe=None)
-    ax1 = fig.add_subplot(gs[0:4,0:4], projection=myproj)
-    ax2 = fig.add_subplot(gs[0:4,5:], projection=myproj)
-    ax3 = fig.add_subplot(gs[5:,1:8], projection=myproj)
-    ax1.set_extent(extent)
-    ax2.set_extent(extent)
-    ax3.set_extent(extent)
-
-    fline_wd = 0.5  # line width
-    falpha = 0.3    # transparency
-
-  # natural_earth
-#    land=cfeature.NaturalEarthFeature('physical','land',back_res,
-#                      edgecolor='face',facecolor=cfeature.COLORS['land'],
-#                      alpha=falpha)
-    lakes=cfeature.NaturalEarthFeature('physical','lakes',back_res,
-                      edgecolor='blue',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-    coastline=cfeature.NaturalEarthFeature('physical','coastline',
-                      back_res,edgecolor='blue',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-    states=cfeature.NaturalEarthFeature('cultural','admin_1_states_provinces',
-                      back_res,edgecolor='black',facecolor='none',
-                      linewidth=fline_wd,linestyle=':',alpha=falpha)
-    borders=cfeature.NaturalEarthFeature('cultural','admin_0_countries',
-                      back_res,edgecolor='red',facecolor='none',
-                      linewidth=fline_wd,alpha=falpha)
-
-  # All lat lons are earth relative, so setup the associated projection correct for that data
-    transform = ccrs.PlateCarree()
-
-  # high-resolution background images
-    if back_img=='on':
-       img = plt.imread(CARTOPY_DIR+'/raster_files/NE1_50M_SR_W.tif')
-       ax1.imshow(img, origin='upper', transform=transform)
-       ax2.imshow(img, origin='upper', transform=transform)
-       ax3.imshow(img, origin='upper', transform=transform)
-
-#  ax.add_feature(land)
-    ax1.add_feature(lakes)
-    ax1.add_feature(states)
-    ax1.add_feature(borders)
-    ax1.add_feature(coastline)
-    ax2.add_feature(lakes)
-    ax2.add_feature(states)
-    ax2.add_feature(borders)
-    ax2.add_feature(coastline)
-    ax3.add_feature(lakes)
-    ax3.add_feature(states)
-    ax3.add_feature(borders)
-    ax3.add_feature(coastline)
-
-# Map/figure has been set up here, save axes instances for use again later
-    keep_ax_lst_1 = ax1.get_children()[:]
-    keep_ax_lst_2 = ax2.get_children()[:]
-    keep_ax_lst_3 = ax3.get_children()[:]
-
-# colors for difference plots, only need to define once
-    diffcolors = ['blue','#1874CD','dodgerblue','deepskyblue','turquoise','white','white','#EEEE00','#EEC900','darkorange','orangered','red']
-
-
-################################
-  # Plot SLP
-################################
-    t1 = time.perf_counter()
-    print(('Working on slp for '+dom))
-
-    units = 'mb'
-    clevs = [976,980,984,988,992,996,1000,1004,1008,1012,1016,1020,1024,1028,1032,1036,1040,1044,1048,1052]
-    clevsdiff = [-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12]
-    cm = plt.cm.Spectral_r
-    cmdiff = matplotlib.colors.ListedColormap(diffcolors)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs1_a = ax1.pcolormesh(lon_shift,lat_shift,slp_1,transform=transform,cmap=cm,norm=norm)
-    cbar1 = plt.colorbar(cs1_a,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=5)
-    cs1_b = ax1.contour(lon_shift,lat_shift,slpsmooth_1,np.arange(940,1060,4),colors='black',linewidths=1.25,transform=transform)
-    plt.clabel(cs1_b,np.arange(940,1060,4),inline=1,fmt='%d',fontsize=6)
-    ax1.text(.5,1.03,'FV3-LAM SLP ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs2_a = ax2.pcolormesh(lon2_shift,lat2_shift,slp_2,transform=transform,cmap=cm,norm=norm)
-    cbar2 = plt.colorbar(cs2_a,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=5)
-    cs2_b = ax2.contour(lon2_shift,lat2_shift,slpsmooth_2,np.arange(940,1060,4),colors='black',linewidths=1.25,transform=transform)
-    plt.clabel(cs2_b,np.arange(940,1060,4),inline=1,fmt='%d',fontsize=6)
-    ax2.text(.5,1.03,'FV3-LAM-2 SLP ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,slp_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=5)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM SLP ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/slp_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot slp for: '+dom) % t3)
-
-
-#################################
-  # Plot 2-m T
-#################################
-    t1 = time.perf_counter()
-    print(('Working on t2m for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = '\xb0''F'
-    clevs = np.linspace(-16,134,51)
-    clevsdiff = [-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6]
-#  cm = plt.cm.Spectral_r
-    cm = cmap_t2m()
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,tmp2m_1,transform=transform,cmap=cm,norm=norm)
-    cs_1.cmap.set_under('white')
-    cs_1.cmap.set_over('white')
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=[-16,-4,8,20,32,44,56,68,80,92,104,116,128],extend='both')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=5)
-    ax1.text(.5,1.03,'FV3-LAM 2-m Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,tmp2m_2,transform=transform,cmap=cm,norm=norm)
-    cs_2.cmap.set_under('white')
-    cs_2.cmap.set_over('white')
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=[-16,-4,8,20,32,44,56,68,80,92,104,116,128],extend='both')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=5)
-    ax2.text(.5,1.03,'FV3-LAM-2 2-m Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,tmp2m_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 2-m Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/2mt_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 2mt for: '+dom) % t3)
-
-
-#################################
-  # Plot 2-m Dew Point
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 2mdew for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = '\xb0''F'
-    clevs = np.linspace(-5,80,35)
-    clevsdiff = [-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12]
-    cm = cmap_q2m()
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,dew2m_1,transform=transform,cmap=cm,norm=norm)
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=6)
-    ax1.text(.5,1.03,'FV3-LAM 2-m Dew Point Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,dew2m_2,transform=transform,cmap=cm,norm=norm)
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=6)
-    ax2.text(.5,1.03,'FV3-LAM-2 2-m Dew Point Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,dew2m_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 2-m Dew Point Temperature ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/2mdew_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 2mdew for: '+dom) % t3)
-
-
-#################################
-  # Plot 10-m WSPD
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 10mwspd for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = 'kts'
-  # Places a wind barb every ~180 km, optimized for CONUS domain
-    skip = round(177.28*(dx/1000.)**-.97)
-    print('skipping every '+str(skip)+' grid points to plot')
-    barblength = 4
-
-    clevs = [5,10,15,20,25,30,35,40,45,50,55,60]
-    clevsdiff = [-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12]
-    colorlist = ['turquoise','dodgerblue','blue','#FFF68F','#E3CF57','peru','brown','crimson','red','fuchsia','DarkViolet']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,wspd10m_1,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=6)
-    ax1.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],uwind_1[::skip,::skip],vwind_1[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax1.text(.5,1.03,'FV3-LAM 10-m Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,wspd10m_2,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_2.cmap.set_under('white',alpha=0.)
-    cs_2.cmap.set_over('black')
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=6)
-    ax2.barbs(lon2_shift[::skip,::skip],lat2_shift[::skip,::skip],uwind_2[::skip,::skip],vwind_2[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax2.text(.5,1.03,'FV3-LAM-2 10-m Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,wspd10m_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 10-m Winds ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/10mwind_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 10mwspd for: '+dom) % t3)
-
-
-#################################
-  # Plot Surface-Based CAPE/CIN
-#################################
-    t1 = time.perf_counter()
-    print(('Working on surface-based CAPE/CIN for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = 'J/kg'
-    clevs = [100,250,500,1000,1500,2000,2500,3000,3500,4000,4500,5000]
-    clevs2 = [-2000,-500,-250,-100,-25]
-    clevsdiff = [-2000,-1500,-1000,-500,-250,-100,0,100,250,500,1000,1500,2000]
-    colorlist = ['blue','dodgerblue','cyan','mediumspringgreen','#FAFAD2','#EEEE00','#EEC900','darkorange','crimson','darkred','darkviolet']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,cape_1,transform=transform,cmap=cm,vmin=100,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=4)
-    cs_1b = ax1.contourf(lon_shift,lat_shift,cin_1,clevs2,colors='none',hatches=['**','++','////','..'],transform=transform)
-    ax1.text(.5,1.05,'FV3-LAM Surface-Based CAPE (shaded) and CIN (hatched) ('+units+') \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,cape_2,transform=transform,cmap=cm,vmin=100,norm=norm)
-    cs_2.cmap.set_under('white',alpha=0.)
-    cs_2.cmap.set_over('black')
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=4)
-    cs_2b = ax2.contourf(lon2_shift,lat2_shift,cin_2,clevs2,colors='none',hatches=['**','++','////','..'],transform=transform)
-    ax2.text(.5,1.05,'FV3-LAM-2 Surface-Based CAPE (shaded) and CIN (hatched) ('+units+') \n <-500 (*), -500<-250 (+), -250<-100 (/), -100<-25 (.) \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,cape_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM Surface-Based CAPE ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/sfcape_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot surface-based CAPE/CIN for: '+dom) % t3)
-
-
-#################################
-  # Plot 500 mb HGT/WIND/VORT
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 500 mb Hgt/Wind/Vort for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = 'x10${^5}$ s${^{-1}}$'
-    skip = round(177.28*(dx/1000.)**-.97)
-    print('skipping every '+str(skip)+' grid points to plot')
-    barblength = 4
-
-    vortlevs = [16,20,24,28,32,36,40]
-    clevsdiff = [-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6]
-    colorlist = ['yellow','gold','goldenrod','orange','orangered','red']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(vortlevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs1_a = ax1.pcolormesh(lon_shift,lat_shift,vort500_1,transform=transform,cmap=cm,norm=norm)
-    cs1_a.cmap.set_under('white')
-    cs1_a.cmap.set_over('darkred')
-    cbar1 = plt.colorbar(cs1_a,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=vortlevs,extend='both')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=6)
-    ax1.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],u500_1[::skip,::skip],v500_1[::skip,::skip],length=barblength,linewidth=0.5,color='steelblue',transform=transform)
-    cs1_b = ax1.contour(lon_shift,lat_shift,z500_1,np.arange(486,600,6),colors='black',linewidths=1,transform=transform)
-    plt.clabel(cs1_b,np.arange(486,600,6),inline_spacing=1,fmt='%d',fontsize=8)
-    ax1.text(.5,1.03,'FV3-LAM 500 mb Heights (dam), Winds (kts), and $\zeta$ ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs2_a = ax2.pcolormesh(lon2_shift,lat2_shift,vort500_2,transform=transform,cmap=cm,norm=norm)
-    cs2_a.cmap.set_under('white')
-    cs2_a.cmap.set_over('darkred')
-    cbar2 = plt.colorbar(cs2_a,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=vortlevs,extend='both')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=6)
-    ax2.barbs(lon2_shift[::skip,::skip],lat2_shift[::skip,::skip],u500_2[::skip,::skip],v500_2[::skip,::skip],length=barblength,linewidth=0.5,color='steelblue',transform=transform)
-    cs2_b = ax2.contour(lon2_shift,lat2_shift,z500_2,np.arange(486,600,6),colors='black',linewidths=1,transform=transform)
-    plt.clabel(cs2_b,np.arange(486,600,6),inline_spacing=1,fmt='%d',fontsize=8)
-    ax2.text(.5,1.03,'FV3-LAM-2 500 mb Heights (dam), Winds (kts), and $\zeta$ ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=8,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,z500_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 500-mb Heights (dam) \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/500_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 500 mb Hgt/Wind/Vort for: '+dom) % t3)
-
-
-#################################
-  # Plot 250 mb WIND
-#################################
-    t1 = time.perf_counter()
-    print(('Working on 250 mb WIND for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = 'kts'
-    skip = round(177.28*(dx/1000.)**-.97)
-    print('skipping every '+str(skip)+' grid points to plot')
-    barblength = 4
-
-    clevs = [50,60,70,80,90,100,110,120,130,140,150]
-    clevsdiff = [-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30]
-    colorlist = ['turquoise','deepskyblue','dodgerblue','#1874CD','blue','beige','khaki','peru','brown','crimson']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-    normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,wspd250_1,transform=transform,cmap=cm,vmin=50,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('red')
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=6)
-    ax1.barbs(lon_shift[::skip,::skip],lat_shift[::skip,::skip],u250_1[::skip,::skip],v250_1[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax1.text(.5,1.03,'FV3-LAM 250 mb Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,wspd250_2,transform=transform,cmap=cm,vmin=50,norm=norm)
-    cs_2.cmap.set_under('white',alpha=0.)
-    cs_2.cmap.set_over('red')
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=6)
-    ax2.barbs(lon2_shift[::skip,::skip],lat2_shift[::skip,::skip],u250_2[::skip,::skip],v250_2[::skip,::skip],length=barblength,linewidth=0.5,color='black',transform=transform)
-    ax2.text(.5,1.03,'FV3-LAM-2 250 mb Winds ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs = ax3.pcolormesh(lon2_shift,lat2_shift,wspd250_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-    cs.cmap.set_under('darkblue')
-    cs.cmap.set_over('darkred')
-    cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-    cbar3.set_label(units,fontsize=6)
-    cbar3.ax.tick_params(labelsize=6)
-    ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 250-mb Winds ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/250wind_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot 250 mb WIND for: '+dom) % t3)
-
-
-#################################
-  # Plot Total QPF
-#################################
-    if (fhr > 0):		# Do not make total QPF plot for forecast hour 0
-      t1 = time.perf_counter()
-      print(('Working on total qpf for '+dom))
-
-    # Clear off old plottables but keep all the map info
-      cbar1.remove()
-      cbar2.remove()
-      cbar3.remove()
-      clear_plotables(ax1,keep_ax_lst_1,fig)
-      clear_plotables(ax2,keep_ax_lst_2,fig)
-      clear_plotables(ax3,keep_ax_lst_3,fig)
-
-      units = 'in'
-      clevs = [0.01,0.1,0.25,0.5,0.75,1,1.25,1.5,1.75,2,2.5,3,4,5,7,10,15,20]
-      clevsdiff = [-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3]
-      colorlist = ['chartreuse','limegreen','green','blue','dodgerblue','deepskyblue','cyan','mediumpurple','mediumorchid','darkmagenta','darkred','crimson','orangered','darkorange','goldenrod','gold','yellow']
-      cm = matplotlib.colors.ListedColormap(colorlist)
-      norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-      normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-      cs_1 = ax1.pcolormesh(lon_shift,lat_shift,qpf_1,transform=transform,cmap=cm,vmin=0.01,norm=norm)
-      cs_1.cmap.set_under('white',alpha=0.)
-      cs_1.cmap.set_over('pink')
-      cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,extend='max')
-      cbar1.set_label(units,fontsize=6)
-      cbar1.ax.set_xticklabels([0.1,0.5,1,1.5,2,3,5,10,20])
-      cbar1.ax.tick_params(labelsize=6)
-      ax1.text(.5,1.03,'FV3-LAM '+fhour+'-hr Accumulated Precipitation ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,qpf_2,transform=transform,cmap=cm,vmin=0.01,norm=norm)
-      cs_2.cmap.set_under('white',alpha=0.)
-      cs_2.cmap.set_over('pink')
-      cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,extend='max')
-      cbar2.set_label(units,fontsize=6)
-      cbar2.ax.set_xticklabels([0.1,0.5,1,1.5,2,3,5,10,20])
-      cbar2.ax.tick_params(labelsize=6)
-      ax2.text(.5,1.03,'FV3-LAM-2 '+fhour+'-hr Accumulated Precipitation ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      cs = ax3.pcolormesh(lon2_shift,lat2_shift,qpf_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-      cs.cmap.set_under('darkblue')
-      cs.cmap.set_over('darkred')
-      cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-      cbar3.set_label(units,fontsize=6)
-      cbar3.ax.tick_params(labelsize=6)
-      ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM '+fhour+'-hr Accumulated Precipitation ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/qpf_diff_'+dom+'_f'+fhour+'.png')
-      t2 = time.perf_counter()
-      t3 = round(t2-t1, 3)
-      print(('%.3f seconds to plot total qpf for: '+dom) % t3)
-
-
-#################################
-  # Plot Max/Min Hourly 2-5 km UH
-#################################
-# Do not make max/min hourly 2-5 km UH plot for forecast hour 0
-      t1 = time.perf_counter()
-      print(('Working on Max/Min Hourly 2-5 km UH for '+dom))
-
-    # Clear off old plottables but keep all the map info
-      cbar1.remove()
-      cbar2.remove()
-      cbar3.remove()
-      clear_plotables(ax1,keep_ax_lst_1,fig)
-      clear_plotables(ax2,keep_ax_lst_2,fig)
-      clear_plotables(ax3,keep_ax_lst_3,fig)
-
-      units = 'm${^2}$ s$^{-2}$'
-      clevs = [-150,-100,-75,-50,-25,-10,0,10,25,50,75,100,150,200,250,300]
-#   alternative colormap for just max UH if you don't want to plot the min UH too
-#   colorlist = ['white','skyblue','mediumblue','green','orchid','firebrick','#EEC900','DarkViolet']
-      colorlist = ['blue','#1874CD','dodgerblue','deepskyblue','turquoise','#E5E5E5','#E5E5E5','#EEEE00','#EEC900','darkorange','orangered','red','firebrick','mediumvioletred','darkviolet']
-      cm = matplotlib.colors.ListedColormap(colorlist)
-      cmdiff = matplotlib.colors.ListedColormap(diffcolors)
-      norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-      normdiff = matplotlib.colors.BoundaryNorm(clevsdiff, cmdiff.N)
-
-      cs_1 = ax1.pcolormesh(lon_shift,lat_shift,uh25_1,transform=transform,cmap=cm,norm=norm)
-      cs_1.cmap.set_under('darkblue')
-      cs_1.cmap.set_over('black')
-      cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-      cbar1.set_label(units,fontsize=6)
-      cbar1.ax.tick_params(labelsize=6)
-      ax1.text(.5,1.03,'FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,uh25_2,transform=transform,cmap=cm,norm=norm)
-      cs_2.cmap.set_under('darkblue')
-      cs_2.cmap.set_over('black')
-      cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-      cbar2.set_label(units,fontsize=6)
-      cbar2.ax.tick_params(labelsize=6)
-      ax2.text(.5,1.03,'FV3-LAM-2 1-h Max/Min 2-5 km Updraft Helicity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      cs = ax3.pcolormesh(lon2_shift,lat2_shift,uh25_diff,transform=transform,cmap=cmdiff,norm=normdiff)
-      cs.cmap.set_under('darkblue')
-      cs.cmap.set_over('darkred')
-      cbar3 = plt.colorbar(cs,ax=ax3,orientation='horizontal',pad=0.05,shrink=0.6,extend='both')
-      cbar3.set_label(units,fontsize=6)
-      cbar3.ax.tick_params(labelsize=6)
-      ax3.text(.5,1.03,'FV3-LAM-2 - FV3-LAM 1-h Max/Min 2-5 km Updraft Helicity ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-      compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/uh25_diff_'+dom+'_f'+fhour+'.png')
-      t2 = time.perf_counter()
-      t3 = round(t2-t1, 3)
-      print(('%.3f seconds to plot Max/Min Hourly 2-5 km UH for: '+dom) % t3)
-
-
-#################################
-  # Plot composite reflectivity
-#################################
-    t1 = time.perf_counter()
-    print(('Working on composite reflectivity for '+dom))
-
-  # Clear off old plottables but keep all the map info
-    cbar1.remove()
-    cbar2.remove()
-    cbar3.remove()
-    clear_plotables(ax1,keep_ax_lst_1,fig)
-    clear_plotables(ax2,keep_ax_lst_2,fig)
-    clear_plotables(ax3,keep_ax_lst_3,fig)
-
-    units = 'dBZ'
-    clevs = np.linspace(5,70,14)
-    clevsdiff = [20,1000]
-    colorlist = ['turquoise','dodgerblue','mediumblue','lime','limegreen','green','#EEEE00','#EEC900','darkorange','red','firebrick','darkred','fuchsia']
-    cm = matplotlib.colors.ListedColormap(colorlist)
-    norm = matplotlib.colors.BoundaryNorm(clevs, cm.N)
-
-    cs_1 = ax1.pcolormesh(lon_shift,lat_shift,refc_1,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_1.cmap.set_under('white',alpha=0.)
-    cs_1.cmap.set_over('black')
-    cbar1 = plt.colorbar(cs_1,ax=ax1,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar1.set_label(units,fontsize=6)
-    cbar1.ax.tick_params(labelsize=6)
-    ax1.text(.5,1.03,'FV3-LAM Composite Reflectivity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax1.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    cs_2 = ax2.pcolormesh(lon2_shift,lat2_shift,refc_2,transform=transform,cmap=cm,vmin=5,norm=norm)
-    cs_2.cmap.set_under('white',alpha=0.)
-    cs_2.cmap.set_over('black')
-    cbar2 = plt.colorbar(cs_2,ax=ax2,orientation='horizontal',pad=0.05,shrink=0.6,ticks=clevs,extend='max')
-    cbar2.set_label(units,fontsize=6)
-    cbar2.ax.tick_params(labelsize=6)
-    ax2.text(.5,1.03,'FV3-LAM-2 Composite Reflectivity ('+units+') \n initialized: '+itime+' valid: '+vtime + ' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax2.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    csdiff = ax3.contourf(lon_shift,lat_shift,refc_1,clevsdiff,colors='red',transform=transform)
-    csdiff2 = ax3.contourf(lon2_shift,lat2_shift,refc_2,clevsdiff,colors='dodgerblue',transform=transform)
-    ax3.text(.5,1.03,'FV3-LAM (red) and FV3-LAM-2 (blue) Composite Reflectivity > 20 ('+units+') \n initialized: '+itime+' valid: '+vtime+' (f'+fhour+')',horizontalalignment='center',fontsize=6,transform=ax3.transAxes,bbox=dict(facecolor='white',alpha=0.85,boxstyle='square,pad=0.2'))
-
-    compress_and_save(EXPT_DIR_1+'/'+ymdh+'/postprd/refc_diff_'+dom+'_f'+fhour+'.png')
-    t2 = time.perf_counter()
-    t3 = round(t2-t1, 3)
-    print(('%.3f seconds to plot composite reflectivity for: '+dom) % t3)
-
-
-######################################################
-
-    t3dom = round(t2-t1dom, 3)
-    print(("%.3f seconds to plot all variables for: "+dom) % t3dom)
-    plt.clf()
-
-######################################################
-
-  main()
-
diff --git a/ush/Python/qsub_job.sh b/ush/Python/qsub_job.sh
deleted file mode 100755
index 6a6b0833f1..0000000000
--- a/ush/Python/qsub_job.sh
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/bin/sh
-#PBS -A an_account 
-#PBS -q regular
-#PBS -l select=1:mpiprocs=24:ncpus=24
-#PBS -l walltime=00:20:00
-#PBS -N plot_allvars
-#PBS -j oe -o plot_allvars.out 
-
-# Prior to submitting the script the following environment variables
-# must be set using export or setenv
-# HOMEdir=/path-to/ufs-srweather-app/regional_workflow
-# EXPTDIR=/path-to/expt_dirs/your_experiment
-
-cd ${HOMEdir}/ush/Python
-set -x
-
-source /etc/profile.d/modules.sh
-
-############
-# Python environment for Cheyenne 
-############
-module load ncarenv
-module load conda/latest
-conda activate /glade/p/ral/jntp/UFS_SRW_app/conda/python_graphics
-
-############
-# Path to shape files
-############
-#Cheyenne:
-SHAPE_FILES=/glade/p/ral/jntp/UFS_SRW_app/tools/NaturalEarth
-
-export GLOBAL_VAR_DEFNS_FP="${EXPTDIR}/var_defns.sh"
-source ${GLOBAL_VAR_DEFNS_FP}
-export CDATE=${DATE_FIRST_CYCL}
-export FCST_START=3
-export FCST_END=${FCST_LEN_HRS}
-export FCST_INC=3
-
-# Usage statement:      Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour in HHH format
-#                       3. Ending forecast hour in HHH format
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR: Experiment directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR/YYYYMMDDHH/postprd/
-#                       6. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                       7. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecast and in constructing the names
-#                          of the post output files.
-
-
-python plot_allvars.py ${CDATE} ${FCST_START} ${FCST_END} ${FCST_INC} \
-                       ${EXPTDIR} ${SHAPE_FILES} ${POST_OUTPUT_DOMAIN_NAME}
diff --git a/ush/Python/qsub_job_diff.sh b/ush/Python/qsub_job_diff.sh
deleted file mode 100755
index e024958203..0000000000
--- a/ush/Python/qsub_job_diff.sh
+++ /dev/null
@@ -1,63 +0,0 @@
-#!/bin/sh
-#PBS -A an_account 
-#PBS -q regular
-#PBS -l select=1:mpiprocs=24:ncpus=24
-#PBS -l walltime=00:30:00
-#PBS -N plot_allvars_diff
-#PBS -j oe -o plot_allvars_diff.out 
-
-# Prior to submitting the script the following environment variables
-# must be set using export or setenv
-# HOMEdir=/path-to/ufs-srweather-app/regional_workflow
-# EXPTDIR1=/path-to/expt_dirs/your_experiment1
-# EXPTDIR2=/path-to/expt_dirs/your_experiment2
-
-cd ${HOMEdir}/ush/Python
-set -x
-
-source /etc/profile.d/modules.sh
-
-############
-# Python environment for Cheyenne 
-############
-module load ncarenv
-module load conda/latest
-conda activate /glade/p/ral/jntp/UFS_SRW_app/conda/python_graphics
-
-############
-# Path to shape files
-############
-#Cheyenne:
-SHAPE_FILES=/glade/p/ral/jntp/UFS_SRW_app/tools/NaturalEarth
-
-export GLOBAL_VAR_DEFNS_FP="${EXPTDIR1}/var_defns.sh"
-source ${GLOBAL_VAR_DEFNS_FP}
-export CDATE=${DATE_FIRST_CYCL}
-export FCST_START=3
-export FCST_END=${FCST_LEN_HRS}
-export FCST_INC=3
-
-# Usage statement:      Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour in HHH format
-#                       3. Ending forecast hour in HHH format
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR_1: Experiment 1 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_1/YYYYMMDDHH/postprd/
-#                       6. EXPT_DIR_2: Experiment 2 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_2/YYYYMMDDHH/postprd/
-#                       7. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                       8. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecasts and in constructing the names 
-#                          of the post output files.  This must be the same 
-#                          for both forecasts.
-
-python plot_allvars_diff.py ${CDATE} ${FCST_START} ${FCST_END} ${FCST_INC} \
-                            ${EXPTDIR1} ${EXPTDIR2} \
-                            ${SHAPE_FILES} \
-                            ${POST_OUTPUT_DOMAIN_NAME}
diff --git a/ush/Python/sq_job.sh b/ush/Python/sq_job.sh
deleted file mode 100755
index 0c06369d29..0000000000
--- a/ush/Python/sq_job.sh
+++ /dev/null
@@ -1,78 +0,0 @@
-#!/bin/sh
-
-#SBATCH --account=an_account
-#SBATCH --qos=batch
-#SBATCH --ntasks=4
-#SBATCH --time=0:20:00
-#SBATCH --job-name="plot_allvars"
-#SBATCH --out=plot_allvars.out
-
-# Prior to submitting the script the following environment variables
-# must be set using export or setenv
-# HOMEdir=/path-to/ufs-srweather-app/regional_workflow
-# EXPTDIR=/path-to/expt_dirs/your_experiment
-
-cd ${HOMEdir}/ush/Python
-set -x
-. /apps/lmod/lmod/init/sh
-
-module purge
-module load hpss
-
-############
-# Python environment for Jet and Hera
-module use -a /contrib/miniconda3/modulefiles
-module load miniconda3
-conda activate pygraf
-
-############
-# Python environment for Orion
-############
-#module use -a /apps/contrib/miniconda3-noaa-gsl/modulefiles
-#module load miniconda3
-#conda activate pygraf
-
-############
-# Python environment for Gaea
-############
-#module use /lustre/f2/pdata/esrl/gsd/contrib/modulefiles
-#module load miniconda3/4.8.3-regional-workflow
-
-############
-# Path to shape files
-############
-#Hera:
-SHAPE_FILES=/scratch2/BMC/det/UFS_SRW_app/v1p0/fix_files/NaturalEarth
-#Jet: 
-#SHAPE_FILES=/lfs4/BMC/wrfruc/FV3-LAM/NaturalEarth
-#Orion: 
-#SHAPE_FILES=/work/noaa/gsd-fv3-dev/UFS_SRW_App/v1p0/fix_files/NaturalEarth
-#Gaea: 
-#SHAPE_FILES=/lustre/f2/pdata/esrl/gsd/ufs/NaturalEarth
-
-export GLOBAL_VAR_DEFNS_FP="${EXPTDIR}/var_defns.sh"
-source ${GLOBAL_VAR_DEFNS_FP}
-export CDATE=${DATE_FIRST_CYCL}
-export FCST_START=6
-export FCST_END=${FCST_LEN_HRS}
-export FCST_INC=6
-
-# Usage statement:      Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour in HHH format
-#                       3. Ending forecast hour in HHH format
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR: Experiment directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR/YYYYMMDDHH/postprd/
-#                       6. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                       7. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecast and in constructing the names
-#                          of the post output files.
-
-
-python plot_allvars.py ${CDATE} ${FCST_START} ${FCST_END} ${FCST_INC} \
-                       ${EXPTDIR} ${SHAPE_FILES} ${POST_OUTPUT_DOMAIN_NAME}
diff --git a/ush/Python/sq_job_diff.sh b/ush/Python/sq_job_diff.sh
deleted file mode 100755
index 5f66a3d7d1..0000000000
--- a/ush/Python/sq_job_diff.sh
+++ /dev/null
@@ -1,85 +0,0 @@
-#!/bin/sh
-
-#SBATCH --account=an_account
-#SBATCH --qos=batch
-#SBATCH --ntasks=8
-#SBATCH --time=0:30:00
-#SBATCH --job-name="plot_allvars_diff"
-#SBATCH --out=plot_allvars.out
-
-# Prior to submitting the script the following environment variables
-# must be set using export or setenv
-# HOMEdir=/path-to/ufs-srweather-app/regional_workflow
-# EXPTDIR1=/path-to/expt_dirs/your_experiment1
-# EXPTDIR2=/path-to/expt_dirs/your_experiment2
-
-cd ${HOMEdir}/ush/Python
-set -x
-. /apps/lmod/lmod/init/sh
-
-module purge
-module load hpss
-
-############
-# Python environment for Jet and Hera
-############
-module use -a /contrib/miniconda3/modulefiles
-module load miniconda3
-conda activate pygraf
-
-############
-# Python environment for Orion
-############
-#module use -a /apps/contrib/miniconda3-noaa-gsl/modulefiles
-#module load miniconda3
-#conda activate pygraf
-
-############
-# Python environment for Gaea
-############
-#module use /lustre/f2/pdata/esrl/gsd/contrib/modulefiles
-#module load miniconda3/4.8.3-regional-workflow
-
-############
-# Path to shape files
-############
-#Hera:
-SHAPE_FILES=/scratch2/BMC/det/UFS_SRW_app/v1p0/fix_files/NaturalEarth
-#Jet: 
-#SHAPE_FILES=/lfs4/BMC/wrfruc/FV3-LAM/NaturalEarth
-#Orion: 
-#SHAPE_FILES=/work/noaa/gsd-fv3-dev/UFS_SRW_App/v1p0/fix_files/NaturalEarth
-#Gaea: 
-#SHAPE_FILES=/lustre/f2/pdata/esrl/gsd/ufs/NaturalEarth
-
-export GLOBAL_VAR_DEFNS_FP="${EXPTDIR1}/var_defns.sh"
-source ${GLOBAL_VAR_DEFNS_FP}
-export CDATE=${DATE_FIRST_CYCL}
-export FCST_START=6
-export FCST_END=${FCST_LEN_HRS}
-export FCST_INC=3
-
-# Usage statement:      Make sure all the necessary modules can be imported.
-#                       The following command line arguments are needed:
-#                       1. Cycle date/time in YYYYMMDDHH format
-#                       2. Starting forecast hour
-#                       3. Ending forecast hour
-#                       4. Forecast hour increment
-#                       5. EXPT_DIR_1: Experiment 1 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_1/YYYYMMDDHH/postprd/
-#                       6. EXPT_DIR_2: Experiment 2 directory
-#                          -Postprocessed data should be found in the directory:
-#                            EXPT_DIR_2/YYYYMMDDHH/postprd/
-#                       7. CARTOPY_DIR:  Base directory of cartopy shapefiles
-#                          -File structure should be:
-#                            CARTOPY_DIR/shapefiles/natural_earth/cultural/*.shp
-#                       8. POST_OUTPUT_DOMAIN_NAME:  Name of native domain
-#                          used in forecasts and in constructing the names 
-#                          of the post output files.  This must be the same 
-#                          for both forecasts.
-
-python plot_allvars_diff.py ${CDATE} ${FCST_START} ${FCST_END} ${FCST_INC} \
-                            ${EXPTDIR1} ${EXPTDIR2} \
-                            ${SHAPE_FILES} \
-                            ${POST_OUTPUT_DOMAIN_NAME}
diff --git a/ush/config.community.yaml b/ush/config.community.yaml
index 3769bd3183..f211c74a1f 100644
--- a/ush/config.community.yaml
+++ b/ush/config.community.yaml
@@ -13,6 +13,7 @@ platform:
   MRMS_OBS_DIR: ""
   NDAS_OBS_DIR: ""
 workflow:
+  USE_CRON_TO_RELAUNCH: false
   EXPT_SUBDIR: test_community
   CCPP_PHYS_SUITE: FV3_GFS_v16
   PREDEF_GRID_NAME: RRFS_CONUS_25km
@@ -33,6 +34,7 @@ workflow_switches:
   RUN_TASK_VX_POINTSTAT: false
   RUN_TASK_VX_ENSGRID: false
   RUN_TASK_VX_ENSPOINT: false
+  RUN_TASK_PLOT_ALLVARS: false
 task_get_extrn_ics:
   EXTRN_MDL_NAME_ICS: FV3GFS
   FV3GFS_FILE_FMT_ICS: grib2
@@ -43,6 +45,8 @@ task_get_extrn_lbcs:
 task_run_fcst:
   WTIME_RUN_FCST: 02:00:00
   QUILTING: true
+task_plot_allvars:
+  COMOUT_REF: ""
 global:
   DO_ENSEMBLE: false
   NUM_ENS_MEMBERS: 2
diff --git a/ush/config.nco.yaml b/ush/config.nco.yaml
index f6971ac732..26171f2f3d 100644
--- a/ush/config.nco.yaml
+++ b/ush/config.nco.yaml
@@ -6,6 +6,7 @@ user:
   MACHINE: hera
   ACCOUNT: an_account
 workflow:
+  USE_CRON_TO_RELAUNCH: false
   EXPT_SUBDIR: test_nco
   CCPP_PHYS_SUITE: FV3_GFS_v16
   PREDEF_GRID_NAME: RRFS_CONUS_25km
@@ -19,6 +20,7 @@ workflow_switches:
   RUN_TASK_MAKE_GRID: false
   RUN_TASK_MAKE_OROG: false
   RUN_TASK_MAKE_SFC_CLIMO: false
+  RUN_TASK_PLOT_ALLVARS: false
 nco:
   NET: rrfs
   model_ver: v1.0
@@ -35,5 +37,7 @@ task_run_fcst:
   WTIME_RUN_FCST: 01:00:00
   WRITE_DOPOST: true
   QUILTING: true
+task_plot_allvars:
+  COMOUT_REF: ""
 task_run_post:
   POST_OUTPUT_DOMAIN_NAME: conus_25km
diff --git a/ush/config_defaults.yaml b/ush/config_defaults.yaml
index 295fc17781..b6ab952225 100644
--- a/ush/config_defaults.yaml
+++ b/ush/config_defaults.yaml
@@ -361,6 +361,7 @@ platform:
   FIXlut: ""
   FIXorg: ""
   FIXsfc: ""
+  FIXshp: ""
   #
   #-----------------------------------------------------------------------
   #
@@ -857,6 +858,9 @@ workflow_switches:
   # to be run. If this flag is set, both ensemble-stat point verification
   # and point verification of ensemble-stat output is computed.
   #
+  # RUN_TASK_PLOT_ALLVARS:
+  # Flag that determines whether to run python plotting scripts
+  #
   #-----------------------------------------------------------------------
   #
   RUN_TASK_MAKE_GRID: true
@@ -878,6 +882,7 @@ workflow_switches:
   RUN_TASK_VX_ENSGRID: false
   RUN_TASK_VX_ENSPOINT: false
 
+  RUN_TASK_PLOT_ALLVARS: false
 #----------------------------
 # MAKE GRID config parameters
 #-----------------------------
@@ -2019,6 +2024,34 @@ task_run_vx_enspoint_prob:
   WTIME_VX_ENSPOINT_PROB: 01:00:00
   MAXTRIES_VX_ENSPOINT_PROB: 1
 
+#----------------------------
+# PLOT_ALLVARS config parameters
+#-----------------------------
+task_plot_allvars:
+  PLOT_ALLVARS_TN: "plot_allvars"
+  NNODES_PLOT_ALLVARS: 1
+  PPN_PLOT_ALLVARS: 24
+  WTIME_PLOT_ALLVARS: 01:00:00
+  MAXTRIES_PLOT_ALLVARS: 1
+  #-------------------------------------------------------------------------
+  # Reference experiment's COMOUT directory. This is where the GRIB2 files 
+  # from postprocessing are located. Make this a template to compare
+  # multiple cycle and dates. COMOUT_REF should end with:
+  #    nco mode: $PDY/$cyc
+  #    community mode: $PDY$cyc/postprd
+  # We don't do this inside the code, so that we can compare nco vs com runs.
+  #-------------------------------------------------------------------------
+  COMOUT_REF: ""
+  #------------------------------
+  # Plot fcts start and increment
+  #------------------------------
+  PLOT_FCST_START: 0
+  PLOT_FCST_INC: 3
+  #-----------------------------------
+  # By default the end is FCST_LEN_HRS
+  #-----------------------------------
+  PLOT_FCST_END: ""
+
 
 #----------------------------
 # global config parameters
diff --git a/ush/machine/cheyenne.yaml b/ush/machine/cheyenne.yaml
index 932b464937..ed10c964d8 100644
--- a/ush/machine/cheyenne.yaml
+++ b/ush/machine/cheyenne.yaml
@@ -26,5 +26,7 @@ platform:
   FIXlut: /glade/p/ral/jntp/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /glade/p/ral/jntp/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /glade/p/ral/jntp/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /glade/p/ral/jntp/UFS_SRW_App/develop/NaturalEarth
 data:
-  FV3GFS: /glade/p/ral/jntp/UFS_CAM/COMGFS/gfs.${yyyymmdd}/${hh}
+  ics_lbcs:
+    FV3GFS: /glade/p/ral/jntp/UFS_CAM/COMGFS/gfs.${yyyymmdd}/${hh}
diff --git a/ush/machine/gaea.yaml b/ush/machine/gaea.yaml
index 63a50f8d2d..07c070b978 100644
--- a/ush/machine/gaea.yaml
+++ b/ush/machine/gaea.yaml
@@ -26,5 +26,7 @@ platform:
   FIXlut: /lustre/f2/pdata/ncep_shared/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /lustre/f2/pdata/ncep_shared/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /lustre/f2/pdata/ncep_shared/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /lustre/f2/pdata/ncep_shared/UFS_SRW_App/develop/NaturalEarth
 data:
-  FV3GFS: /lustre/f2/dev/Mark.Potts/EPIC/SRW/model_data/FV3GFS/${yyyymmdd}${hh}
+  ics_lbcs:
+    FV3GFS: /lustre/f2/dev/Mark.Potts/EPIC/SRW/model_data/FV3GFS/${yyyymmdd}${hh}
diff --git a/ush/machine/hera.yaml b/ush/machine/hera.yaml
index 3027b1e4fc..3a7769a557 100644
--- a/ush/machine/hera.yaml
+++ b/ush/machine/hera.yaml
@@ -30,4 +30,5 @@ platform:
   FIXlut: /scratch2/BMC/det/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /scratch2/BMC/det/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /scratch2/BMC/det/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /scratch2/BMC/det/UFS_SRW_App/develop/NaturalEarth
   EXTRN_MDL_DATA_STORES: hpss aws nomads
diff --git a/ush/machine/jet.yaml b/ush/machine/jet.yaml
index 535b05672d..22f945335d 100644
--- a/ush/machine/jet.yaml
+++ b/ush/machine/jet.yaml
@@ -30,11 +30,13 @@ platform:
   FIXlut: /mnt/lfs4/BMC/wrfruc/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /mnt/lfs4/BMC/wrfruc/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /mnt/lfs4/BMC/wrfruc/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /mnt/lfs4/BMC/wrfruc/UFS_SRW_App/develop/NaturalEarth
   EXTRN_MDL_DATA_STORES: hpss aws nomads
 data:
-  FV3GFS:
-    nemsio: /public/data/grids/gfs/nemsio
-    grib2: /public/data/grids/gfs/0p25deg/grib2
-    netcdf: /public/data/grids/gfs/anl/netcdf
-  RAP: /public/data/grids/rap/full/wrfprs/grib2
-  HRRR: /public/data/grids/hrrr/conus/wrfprs/grib2
+  ics_lbcs:
+    FV3GFS:
+      nemsio: /public/data/grids/gfs/nemsio
+      grib2: /public/data/grids/gfs/0p25deg/grib2
+      netcdf: /public/data/grids/gfs/anl/netcdf
+    RAP: /public/data/grids/rap/full/wrfprs/grib2
+    HRRR: /public/data/grids/hrrr/conus/wrfprs/grib2
diff --git a/ush/machine/linux.yaml b/ush/machine/linux.yaml
index c7a36ad6bc..db9749aa73 100644
--- a/ush/machine/linux.yaml
+++ b/ush/machine/linux.yaml
@@ -12,5 +12,7 @@ platform:
   FIXlut: /home/username/DATA/UFS/fix/fix_lut
   FIXorg: /home/username/DATA/UFS/fix/fix_orog
   FIXsfc: /home/username/DATA/UFS/fix/fix_sfc_climo
+  FIXshp: /home/username/DATA/UFS/NaturalEarth
 data:
-  FV3GFS: /home/username/DATA/UFS/FV3GFS
+  ics_lbcs:
+    FV3GFS: /home/username/DATA/UFS/FV3GFS
diff --git a/ush/machine/macos.yaml b/ush/machine/macos.yaml
index 506c136b82..7d16746174 100644
--- a/ush/machine/macos.yaml
+++ b/ush/machine/macos.yaml
@@ -12,5 +12,7 @@ platform:
   FIXlut: /Users/username/DATA/UFS/fix/fix_lut
   FIXorg: /Users/username/DATA/UFS/fix/fix_orog
   FIXsfc: /Users/username/DATA/UFS/fix/fix_sfc_climo
+  FIXshp: /Users/username/DATA/UFS/NaturalEarth
 data:
-  FV3GFS: /Users/username/DATA/UFS/FV3GFS
+  ics_lbcs:
+    FV3GFS: /Users/username/DATA/UFS/FV3GFS
diff --git a/ush/machine/noaacloud.yaml b/ush/machine/noaacloud.yaml
index a5800aa7de..61595c5ab3 100644
--- a/ush/machine/noaacloud.yaml
+++ b/ush/machine/noaacloud.yaml
@@ -23,6 +23,8 @@ platform:
   FIXlut: /contrib/EPIC/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /contrib/EPIC/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /contrib/EPIC/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /contrib/EPIC/UFS_SRW_App/develop/NaturalEarth
   EXTRN_MDL_DATA_STORES: aws nomads
 data:
-  FV3GFS: /contrib/GST/model_data/FV3GFS/${yyyymmdd}${hh}
+  ics_lbcs:
+    FV3GFS: /contrib/GST/model_data/FV3GFS/${yyyymmdd}${hh}
diff --git a/ush/machine/odin.yaml b/ush/machine/odin.yaml
index b47a654793..b92c8245b6 100644
--- a/ush/machine/odin.yaml
+++ b/ush/machine/odin.yaml
@@ -24,9 +24,11 @@ platform:
   FIXlut: /scratch/ywang/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /scratch/ywang/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /scratch/ywang/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /scratch/ywang/UFS_SRW_App/develop/NaturalEarth
 data:
-  GSMGFS: /scratch/ywang/UFS_SRW_App/develop/input_model_data/GFS
-  FV3GFS: /scratch/ywang/UFS_SRW_App/develop/input_model_data/FV3GFS
-  HRRR: /scratch/ywang/UFS_SRW_App/develop/input_model_data/HRRR
-  RAP: /scratch/ywang/UFS_SRW_App/develop/input_model_data/RAP
-  NAM: /scratch/ywang/UFS_SRW_App/develop/input_model_data/NAM
+  ics_lbcs:
+    GSMGFS: /scratch/ywang/UFS_SRW_App/develop/input_model_data/GFS
+    FV3GFS: /scratch/ywang/UFS_SRW_App/develop/input_model_data/FV3GFS
+    HRRR: /scratch/ywang/UFS_SRW_App/develop/input_model_data/HRRR
+    RAP: /scratch/ywang/UFS_SRW_App/develop/input_model_data/RAP
+    NAM: /scratch/ywang/UFS_SRW_App/develop/input_model_data/NAM
diff --git a/ush/machine/orion.yaml b/ush/machine/orion.yaml
index a3ca5ba133..c43709c5b0 100644
--- a/ush/machine/orion.yaml
+++ b/ush/machine/orion.yaml
@@ -28,4 +28,5 @@ platform:
   FIXlut: /work/noaa/fv3-cam/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /work/noaa/fv3-cam/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /work/noaa/fv3-cam/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /work/noaa/fv3-cam/UFS_SRW_App/develop/NaturalEarth
   EXTRN_MDL_DATA_STORES: aws nomads
diff --git a/ush/machine/singularity.yaml b/ush/machine/singularity.yaml
index 43399460cd..efa8d9f997 100644
--- a/ush/machine/singularity.yaml
+++ b/ush/machine/singularity.yaml
@@ -20,4 +20,5 @@ platform:
   FIXlut: /contrib/global/glopara/fix/fix_lut
   FIXorg: /contrib/global/glopara/fix/fix_orog
   FIXsfc: /contrib/global/glopara/fix/fix_sfc_climo
+  FIXshp: /contrib/global/glopara/NaturalEarth
   EXTRN_MDL_DATA_STORES: aws nomads
diff --git a/ush/machine/stampede.yaml b/ush/machine/stampede.yaml
index 2e663bbcd2..6200df86bc 100644
--- a/ush/machine/stampede.yaml
+++ b/ush/machine/stampede.yaml
@@ -24,9 +24,11 @@ platform:
   FIXlut: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/NaturalEarth
 data:
-  GSMGFS: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/GFS
-  FV3GFS: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/FV3GFS
-  HRRR: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/HRRR
-  RAP: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/RAP
-  NAM: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/NAM
+  ics_lbcs:
+    GSMGFS: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/GFS
+    FV3GFS: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/FV3GFS
+    HRRR: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/HRRR
+    RAP: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/RAP
+    NAM: /work2/00315/tg455890/stampede2/UFS_SRW_App/develop/input_model_data/NAM
diff --git a/ush/machine/wcoss2.yaml b/ush/machine/wcoss2.yaml
index c304d8cb08..b514f6c386 100644
--- a/ush/machine/wcoss2.yaml
+++ b/ush/machine/wcoss2.yaml
@@ -25,11 +25,13 @@ platform:
   FIXlut: /lfs/h2/emc/lam/noscrub/UFS_SRW_App/develop/fix/fix_lut
   FIXorg: /lfs/h2/emc/lam/noscrub/UFS_SRW_App/develop/fix/fix_orog
   FIXsfc: /lfs/h2/emc/lam/noscrub/UFS_SRW_App/develop/fix/fix_sfc_climo
+  FIXshp: /lfs/h2/emc/lam/noscrub/UFS_SRW_App/develop/NaturalEarth
   EXTRN_MDL_DATA_STORES: hpss aws nomads
 data:
-  GSMGFS: compath.py ${envir}/gsmgfs/${gsmgfs_ver}/gsmgfs.${PDYext}
-  FV3GFS: compath.py ${envir}/gfs/${gfs_ver}/gfs.${PDYext}/${cycext}/atmos
-  RAP: compath.py ${envir}/rap/${rap_ver}/rap.${PDYext}
-  NAM: compath.py ${envir}/nam/${nam_ver}/nam.${PDYext}
-  HRRR: compath.py ${envir}/hrrr/${hrrr_ver}/hrrr.${PDYext}/conus
+  ics_lbcs:
+    GSMGFS: compath.py ${envir}/gsmgfs/${gsmgfs_ver}/gsmgfs.${PDYext}
+    FV3GFS: compath.py ${envir}/gfs/${gfs_ver}/gfs.${PDYext}/${cycext}/atmos
+    RAP: compath.py ${envir}/rap/${rap_ver}/rap.${PDYext}
+    NAM: compath.py ${envir}/nam/${nam_ver}/nam.${PDYext}
+    HRRR: compath.py ${envir}/hrrr/${hrrr_ver}/hrrr.${PDYext}/conus
 
diff --git a/ush/setup.py b/ush/setup.py
index 572433459f..5bc175b33a 100644
--- a/ush/setup.py
+++ b/ush/setup.py
@@ -158,8 +158,12 @@ def setup():
 
     # ics and lbcs
     def get_location(xcs, fmt):
-        if ("data" in machine_cfg) and (xcs in machine_cfg["data"]):
-            v = machine_cfg["data"][xcs]
+        if (
+             "data" in machine_cfg
+             and ("ics_lbcs" in machine_cfg["data"])
+             and (xcs in machine_cfg["data"]["ics_lbcs"])
+        ):
+            v = machine_cfg["data"]["ics_lbcs"][xcs]
             if not isinstance(v, dict):
                 return v
             else: