diff --git a/docs/Users_Guide/config_options.rst b/docs/Users_Guide/config_options.rst
index 35805c71d3..92470cfda5 100644
--- a/docs/Users_Guide/config_options.rst
+++ b/docs/Users_Guide/config_options.rst
@@ -459,12 +459,20 @@ Regions of parallelized code are:
 
   * :code:`fractional_coverage (data_plane_util.cc)`
 
+    * Called by `gen_ens_prod` to compute NMEP outputs.
+    * Called by `grid_stat` when applying neighborhood verification methods.
+
+  * :code:`ShapeData::conv_filter_circ() (shapedata.cc)`
+
+    * Called by `mode` to apply a convolution smoothing operation when
+      defining objects.
+
 Only the following top-level executables can presently benefit from OpenMP
 parallelization:
 
   * :code:`grid_stat`
-  * :code:`ensemble_stat`
   * :code:`grid_ens_prod`
+  * :code:`mode`
 
 **Thread Binding**    
 
@@ -474,7 +482,7 @@ guarantees that threads remain evenly distributed across the available cores.
 Otherwise, the operating system may migrate threads between cores during a run.
 
 OpenMP provides some environment variables to handle this: :code:`OMP_PLACES`
-and  :code:`OMP_PROC_BIND`.  We anticipate that the effect of setting only
+and :code:`OMP_PROC_BIND`.  We anticipate that the effect of setting only
 :code:`OMP_PROC_BIND=true` would be neutral-to-positive.
 
 However, there are sometimes compiler-specific environment variables. Instead,
diff --git a/src/basic/vx_util/data_plane.cc b/src/basic/vx_util/data_plane.cc
index 2dd2dcaf38..52d2b72f9e 100644
--- a/src/basic/vx_util/data_plane.cc
+++ b/src/basic/vx_util/data_plane.cc
@@ -23,6 +23,7 @@
 using namespace std;
 
 #include <algorithm>
+#include <map>
 
 #include "data_plane.h"
 
@@ -230,57 +231,57 @@ void DataPlane::dump(ostream & out, int depth) const {
 
 void DataPlane::debug_examine(bool show_all_values) const {
 
-   vector<double> values;
-   vector<int> count;
+   // Nothing to print if verbosity level is below 4
+   if(mlog.verbosity_level() < 4) return;
+
+   map<double,int> value_count_map;
    int total_count = 0;
-   
-   for (int x=0; x<Nx; ++x) {
-      for (int y=0; y<Ny; ++y) {
-         double v = get(x,y);
-         if (v <= 0) {
-            continue;
-         }
-         ++total_count;
-         if (show_all_values) {
-            vector<double>::iterator vi;
-            vi = find(values.begin(), values.end(), v);
-            if (vi == values.end()) {
-               values.push_back(v);
-               count.push_back(1);
-            } else {
-               int ii = vi - values.begin();
-               count[ii] = count[ii] + 1;
-            }
+
+   for(int i=0; i<Data.size(); i++) {
+
+      // Count positive values
+      if(Data[i] > 0) total_count++;
+
+      if (show_all_values) {
+
+         // Add a new map entry
+         if(value_count_map.count(Data[i]) == 0) {
+            value_count_map[Data[i]] = 0;
          }
+
+         // Increment count
+         value_count_map[Data[i]] += 1;
       }
    }
-   if (show_all_values) {
-      for (size_t i=0; i<values.size(); ++i) {
-         mlog << Debug(4) << " data value=" << values[i] << " count=" << count[i] << "\n";
+
+   if(show_all_values) {
+      for(auto it  = value_count_map.begin();
+               it != value_count_map.end(); it++) {
+         mlog << Debug(4) << " data value=" << it->first
+              << " count=" << it->second << "\n";
       }
    }
-   mlog << Debug(4) << "Total count = " << total_count << "\n";
+
+   mlog << Debug(4) << "Total count > 0 = " << total_count
+        << " of " << Data.size() << "\n";
+
+   return;
 }
 
 ///////////////////////////////////////////////////////////////////////////////
 
 string DataPlane::sdebug_examine() const{
+   ConcatString cs;
+   int n = 0;
 
-   int total_count = 0;
-   
-   for (int x=0; x<Nx; ++x) {
-      for (int y=0; y<Ny; ++y) {
-         double v = get(x,y);
-         if (v <= 0) {
-            continue;
-         }
-         ++total_count;
-      }
+   // Count positive values
+   for(auto it = Data.begin(); it != Data.end(); it++) {
+      if(*it > 0) n++;
    }
-   char buf[1000];
-   sprintf(buf, "Total count = %d", total_count);
-   string retval = buf;
-   return retval;
+
+   cs << "Total count > 0 = " << n << " of " << (int) Data.size();
+
+   return cs;
 }
 
 ///////////////////////////////////////////////////////////////////////////////
diff --git a/src/basic/vx_util/vx_util.h b/src/basic/vx_util/vx_util.h
index 84638028ff..c50e7d16fc 100644
--- a/src/basic/vx_util/vx_util.h
+++ b/src/basic/vx_util/vx_util.h
@@ -37,6 +37,7 @@
 #include "fix_float.h"
 #include "get_filenames.h"
 #include "grib_constants.h"
+#include "GridTemplate.h"
 #include "int_array.h"
 #include "interp_mthd.h"
 #include "interp_util.h"
diff --git a/src/libcode/vx_shapedata/engine.cc b/src/libcode/vx_shapedata/engine.cc
index d16ba7de6c..0bfab73929 100644
--- a/src/libcode/vx_shapedata/engine.cc
+++ b/src/libcode/vx_shapedata/engine.cc
@@ -6,7 +6,6 @@
 // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
 // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 
-
 ///////////////////////////////////////////////////////////////////////
 
 using namespace std;
@@ -17,6 +16,8 @@ using namespace std;
 #include <stdlib.h>
 #include <string.h>
 #include <cmath>
+#include <set>
+#include <map>
 
 #include "engine.h"
 #include "mode_columns.h"
@@ -591,12 +592,17 @@ void ModeFuzzyEngine::do_fcst_convolution() {
    // Apply a circular convolution to the raw field
    //
 
-   mlog << Debug(4) << " Before fcst convolution: " << fcst_conv->sdebug_examine() << "\n";
-      
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " Before fcst convolution: "
+           << fcst_conv->sdebug_examine() << "\n";
+   }
 
    if(r > 0) fcst_conv->conv_filter_circ(2*r + 1, conf_info.Fcst->vld_thresh);
 
-   mlog << Debug(4) << " After fcst convolution: " << fcst_conv->sdebug_examine() << "\n";
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After fcst convolution: "
+           << fcst_conv->sdebug_examine() << "\n";
+   }
 
    need_fcst_conv       = false;
    need_fcst_thresh     = true;
@@ -625,14 +631,20 @@ void ModeFuzzyEngine::do_obs_convolution() {
    mlog << Debug(3) << "Applying circular convolution of radius "
         << r << " to the observation field.\n";
 
-   mlog << Debug(4) << " Before obs convolution: " << obs_conv->sdebug_examine() << "\n";
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " Before obs convolution: "
+           << obs_conv->sdebug_examine() << "\n";
+   }
 
    //
    // Apply a circular convolution to the raw field
    //
    if(r > 0) obs_conv->conv_filter_circ(2*r + 1, conf_info.Obs->vld_thresh);
 
-   mlog << Debug(4) << " After obs convolution: " << obs_conv->sdebug_examine() << "\n";
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After obs convolution: "
+           << obs_conv->sdebug_examine() << "\n";
+   }
 
    need_obs_conv       = false;
    need_obs_thresh     = true;
@@ -657,13 +669,21 @@ void ModeFuzzyEngine::do_fcst_thresholding() {
 
    *fcst_thresh = *fcst_raw;
    fcst_thresh->threshold(conf_info.Fcst->conv_thresh);
-   mlog << Debug(4) << " After thresholding raw fcst field: " << fcst_thresh->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After thresholding raw fcst field: "
+           << fcst_thresh->sdebug_examine() << "\n";
+   }
 
    //
    // Threshold the convolved field
    //
    fcst_mask->threshold(conf_info.Fcst->conv_thresh);
-   mlog << Debug(4) << " After thresholding convolved fcst field: " << fcst_mask->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After thresholding convolved fcst field: "
+           << fcst_mask->sdebug_examine() << "\n";
+   }
 
    mlog << Debug(3) << "Applying convolution threshold "
         << conf_info.Fcst->conv_thresh.get_str()
@@ -692,13 +712,21 @@ void ModeFuzzyEngine::do_obs_thresholding() {
 
    *obs_thresh = *obs_raw;
    obs_thresh->threshold(conf_info.Obs->conv_thresh);
-   mlog << Debug(4) << " After thresholding raw obs field: " << obs_thresh->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After thresholding raw obs field: "
+           << obs_thresh->sdebug_examine() << "\n";
+   }
 
    //
    // Threshold the convolved field
    //
    obs_mask->threshold(conf_info.Obs->conv_thresh);
-   mlog << Debug(4) << " After thresholding convolved obs field: " << obs_mask->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After thresholding convolved obs field: "
+           << obs_mask->sdebug_examine() << "\n";
+   }
 
    mlog << Debug(3) << "Applying convolution threshold "
         << conf_info.Obs->conv_thresh.get_str()
@@ -733,7 +761,11 @@ void ModeFuzzyEngine::do_fcst_filtering() {
                                 fcst_raw, conf_info.Fcst->conv_thresh, grid,
                                 conf_info.Fcst->var_info->is_precipitation());
 
-      mlog << Debug(4) << " After attribute filtering of fcst field: " << fcst_mask->sdebug_examine() << "\n";
+      if(mlog.verbosity_level() >= 4) {
+         mlog << Debug(4) << " After attribute filtering of fcst field: "
+              << fcst_mask->sdebug_examine() << "\n";
+      }
+
       mlog << Debug(3) << "Applying object attribute filtering"
            << " resulted in " <<  fcst_mask->n_objects()
            << " simple forecast objects.\n";
@@ -769,7 +801,11 @@ void ModeFuzzyEngine::do_obs_filtering() {
                                obs_raw, conf_info.Obs->conv_thresh, grid,
                                conf_info.Obs->var_info->is_precipitation());
 
-      mlog << Debug(4) << " After attribute filtering of obs field: " << obs_mask->sdebug_examine() << "\n";
+      if(mlog.verbosity_level() >= 4) {
+         mlog << Debug(4) << " After attribute filtering of obs field: "
+              << obs_mask->sdebug_examine() << "\n";
+      }
+
       mlog << Debug(3) << "Applying object attribute filtering"
            << " resulted in " <<  obs_mask->n_objects()
            << " simple observation objects.\n";
@@ -797,7 +833,11 @@ void ModeFuzzyEngine::do_fcst_splitting() {
    if(!need_fcst_split) return;
 
    *fcst_split = split(*fcst_mask, n_fcst);
-   mlog << Debug(4) << " After splitting of fcst field: " << fcst_split->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After splitting of fcst field: "
+           << fcst_split->sdebug_examine() << "\n";
+   }
 
    need_fcst_split      = false;
    need_fcst_merge      = true;
@@ -816,7 +856,11 @@ void ModeFuzzyEngine::do_obs_splitting() {
    if(!need_obs_split) return;
 
    *obs_split = split(*obs_mask, n_obs);
-   mlog << Debug(4) << " After splitting of obs field: " << obs_split->sdebug_examine() << "\n";
+
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After splitting of obs field: "
+           << obs_split->sdebug_examine() << "\n";
+   }
 
    need_obs_split      = false;
    need_obs_merge      = true;
@@ -1033,32 +1077,22 @@ void ModeFuzzyEngine::do_matching() {
 
 void ModeFuzzyEngine::do_no_match() {
    int j, k, n;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * obs_shape = (ShapeData *) 0;
+   ShapeData cur_shape;
 
    do_fcst_splitting();
    do_obs_splitting();
 
    clear_colors();
 
-   fcst_shape = new ShapeData [n_fcst];
-   obs_shape  = new ShapeData [n_obs];
-
-   if(!fcst_shape || !obs_shape) {
-
-      mlog << Error << "\nModeFuzzyEngine::do_no_match() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
-
    //
    // Do the single features
    //
    fcst_single.set_size(n_fcst);
 
    for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-      fcst_single[j].set(*fcst_raw, *fcst_thresh, fcst_shape[j],
+      cur_shape = select(*fcst_split, j+1);
+      fcst_single[j].set(*fcst_raw, *fcst_thresh,
+                         *fcst_split, cur_shape,
                          conf_info.inten_perc_value,
                          conf_info.Fcst->var_info->is_precipitation());
       fcst_single[j].object_number = j+1;
@@ -1067,8 +1101,9 @@ void ModeFuzzyEngine::do_no_match() {
    obs_single.set_size(n_obs);
 
    for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-      obs_single[j].set(*obs_raw, *obs_thresh, obs_shape[j],
+      cur_shape = select(*obs_split, j+1);
+      obs_single[j].set(*obs_raw, *obs_thresh,
+                        *obs_split, cur_shape,
                         conf_info.inten_perc_value,
                         conf_info.Obs->var_info->is_precipitation());
       obs_single[j].object_number = j+1;
@@ -1112,9 +1147,6 @@ void ModeFuzzyEngine::do_no_match() {
    // Done
    //
 
-   delete [] fcst_shape; fcst_shape = (ShapeData *) 0;
-   delete [] obs_shape;  obs_shape  = (ShapeData *) 0;
-
    return;
 }
 
@@ -1123,32 +1155,22 @@ void ModeFuzzyEngine::do_no_match() {
 void ModeFuzzyEngine::do_match_merge() {
    int j, k, n;
    InterestInfo junkinfo;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * obs_shape = (ShapeData *) 0;
+   ShapeData cur_shape;
 
    do_fcst_splitting();
    do_obs_splitting();
 
    clear_colors();
 
-   fcst_shape = new ShapeData [n_fcst];
-   obs_shape = new ShapeData [n_obs];
-
-   if(!fcst_shape || !obs_shape) {
-
-      mlog << Error << "\nModeFuzzyEngine::do_match_merge() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
-
    //
    // Do the single features
    //
    fcst_single.set_size(n_fcst);
 
    for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-      fcst_single[j].set(*fcst_raw, *fcst_thresh, fcst_shape[j],
+      cur_shape = select(*fcst_split, j+1);
+      fcst_single[j].set(*fcst_raw, *fcst_thresh,
+                         *fcst_split, cur_shape,
                          conf_info.inten_perc_value,
                          conf_info.Fcst->var_info->is_precipitation());
       fcst_single[j].object_number = j+1;
@@ -1157,8 +1179,9 @@ void ModeFuzzyEngine::do_match_merge() {
    obs_single.set_size(n_obs);
 
    for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-      obs_single[j].set(*obs_raw, *obs_thresh, obs_shape[j],
+      cur_shape = select(*obs_split, j+1);
+      obs_single[j].set(*obs_raw, *obs_thresh,
+                        *obs_split, cur_shape,
                         conf_info.inten_perc_value,
                         conf_info.Obs->var_info->is_precipitation());
       obs_single[j].object_number = j+1;
@@ -1277,9 +1300,6 @@ void ModeFuzzyEngine::do_match_merge() {
    // Done
    //
 
-   delete [] fcst_shape; fcst_shape = (ShapeData *) 0;
-   delete [] obs_shape;  obs_shape = (ShapeData *) 0;
-
    return;
 }
 
@@ -1292,12 +1312,10 @@ void ModeFuzzyEngine::do_match_merge() {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_fcst_merge_thresh() {
-   int j, k, x, y;
-   int n_fcst_merge, intersection;
-   int count, first_k;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_fcst_merge;
    ShapeData fcst_merge_mask, fcst_merge_split;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * fcst_merge_shape = (ShapeData *) 0;
 
    do_fcst_splitting();
 
@@ -1318,90 +1336,60 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() {
    fcst_merge_split = split(fcst_merge_mask, n_fcst_merge);
 
    //
-   // Allocate space for all of the simple forecast shapes and
-   // forecast merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   fcst_shape       = new ShapeData [n_fcst];
-   fcst_merge_shape = new ShapeData [n_fcst_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!fcst_shape || !fcst_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_fcst;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<fcst_split->data.nxy(); j++) {
+
+      oid = fcst_split->data.buf()[j];
+      mid = fcst_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple forecast shapes and
-   // forecast merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-   }
-
-   for(j=0; j<n_fcst_merge; j++) {
-      fcst_merge_shape[j] = select(fcst_merge_split, j+1);
+   for(oid=1; oid<=n_fcst; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No forecast threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object! "
+              << "Recommend modifying the \"fcst.merge_thresh\" setting.\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_fcst_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_fcst; k++) {
-
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
-
-         for(x=0; x<fcst_merge_mask.data.nx(); x++) {
-            for(y=0; y<fcst_merge_mask.data.ny(); y++) {
+   for(mid=1; mid<=n_fcst_merge; mid++) {
 
-               if(fcst_merge_shape[j].s_is_on(x, y) &&
-                  fcst_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the fcst object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= fcst area
-         //
-
-         if(intersection >= fcst_shape[k].area()) {
-
-            collection.make_room();
-
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(first_k, -1);
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(*it, -1);
+      }
 
-            count++;
-         }
-      }  // end for k
+      // Increment set counter
+      collection.n_sets++;
 
-      if(count > 0) collection.n_sets++;
-   }  // end for j
-
-   //
-   // Done
-   //
-
-   delete [] fcst_shape;       fcst_shape = (ShapeData *) 0;
-   delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
    return;
 }
@@ -1415,12 +1403,10 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_obs_merge_thresh() {
-   int j, k, x, y;
-   int n_obs_merge, intersection;
-   int count, first_k;
+   const char *method_name = "ModeFuzzyEngine::do_obs_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_obs_merge;
    ShapeData obs_merge_mask, obs_merge_split;
-   ShapeData * obs_shape = (ShapeData *) 0;
-   ShapeData * obs_merge_shape = (ShapeData *) 0;
 
    do_obs_splitting();
 
@@ -1441,328 +1427,243 @@ void ModeFuzzyEngine::do_obs_merge_thresh() {
    obs_merge_split = split(obs_merge_mask, n_obs_merge);
 
    //
-   // Allocate space for all of the simple observation shapes and
-   // observation merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   obs_shape       = new ShapeData [n_obs];
-   obs_merge_shape = new ShapeData [n_obs_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!obs_shape || !obs_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_obs;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<obs_split->data.nxy(); j++) {
+
+      oid = obs_split->data.buf()[j];
+      mid = obs_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple observation shapes and
-   // simple merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-   }
-
-   for(j=0; j<n_obs_merge; j++) {
-      obs_merge_shape[j] = select(obs_merge_split, j+1);
+   for(oid=1; oid<=n_obs; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No observation threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object! "
+              << "Recommend modifying the \"obs.merge_thresh\" setting.\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_obs_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_obs; k++) {
-
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
-
-         for(x=0; x<obs_merge_mask.data.nx(); x++) {
-            for(y=0; y<obs_merge_mask.data.ny(); y++) {
-
-               if(obs_merge_shape[j].s_is_on(x, y) &&
-                  obs_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the obs object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= obs area
-         //
-         if(intersection >= obs_shape[k].area()) {
-
-            collection.make_room();
+   for(mid=1; mid<=n_obs_merge; mid++) {
 
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(-1, first_k);
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-
-            count++;
-         }
-      }  // end for k
-
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(-1, *it);
+      }
 
-   //
-   // Done
-   //
+      // Increment set counter
+      collection.n_sets++;
 
-   delete [] obs_shape;       obs_shape = (ShapeData *) 0;
-   delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
    return;
 }
 
+///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_fcst_merge_thresh(const ShapeData &merge_data) {
-   int j, k, x, y;
-   int n_fcst_merge, intersection;
-   int count, first_k;
-   //ShapeData fcst_merge_mask;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_fcst_merge;
    ShapeData fcst_merge_split;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * fcst_merge_shape = (ShapeData *) 0;
 
    do_fcst_splitting();
 
    //
-   // Define the forecast merge field by applying the specified threshold
-   // to the convolved field
-   //
-   fcst_merge_split = merge_data;
-   n_fcst_merge = fcst_merge_split.n_objects();
-      
-   // fcst_merge_mask = *fcst_conv;
-
-   // //
-   // // Threshold the forecast merge field
-   // //
-   // fcst_merge_mask.threshold(conf_info.Fcst->merge_thresh);
-
-   // //
-   // // Split up the forecast merge field
-   // //
-   // fcst_merge_split = split(fcst_merge_mask, n_fcst_merge);
-
-   //
-   // Allocate space for all of the simple forecast shapes and
-   // forecast merge shapes
+   // Check dimensions
    //
-   fcst_shape       = new ShapeData [n_fcst];
-   fcst_merge_shape = new ShapeData [n_fcst_merge];
-
-   if(!fcst_shape || !fcst_merge_shape) {
-
-      mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> "
-           << "memory allocation error\n\n";
+   if(fcst_split->data.nx() != merge_data.data.nx() ||
+      fcst_split->data.ny() != merge_data.data.ny()) {
+      mlog << Error << "\n" << method_name
+           << "grid dimensions do not match (" << fcst_split->data.nx() << ", "
+           << fcst_split->data.ny() << ") != (" << merge_data.data.nx() << ", "
+           << merge_data.data.ny() << ")!\n\n";
       exit(1);
    }
 
    //
-   // Select all of the simple forecast shapes and
-   // forecast merge shapes
+   // Assume that the input merge data has already been split
    //
-   for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-   }
-
-   for(j=0; j<n_fcst_merge; j++) {
-      fcst_merge_shape[j] = select(fcst_merge_split, j+1);
-   }
+   fcst_merge_split = merge_data;
+   n_fcst_merge = fcst_merge_split.n_objects();
 
    //
-   // Calculate the composite object sets
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   for(j=0; j<n_fcst_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_fcst; k++) {
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
+   // Add map entries for each object
+   for(j=1; j<=n_fcst;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-         for(x=0; x<fcst_merge_split.data.nx(); x++) {
-            for(y=0; y<fcst_merge_split.data.ny(); y++) {
+   // Loop over all the data points
+   for(j=0; j<fcst_split->data.nxy(); j++) {
 
-               if(fcst_merge_shape[j].s_is_on(x, y) &&
-                  fcst_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
+      oid = fcst_split->data.buf()[j];
+      mid = fcst_merge_split.data.buf()[j];
 
-         //
-         // Add to set collection only if the fcst object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= fcst area
-         //
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
 
-         if(intersection >= fcst_shape[k].area()) {
+   } // end for j
 
-            collection.make_room();
+   //
+   // Each simple object should be fully contained in a merge object
+   //
+   for(oid=1; oid<=n_fcst; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No forecast threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object!\n\n";
+         return;
+      }
+   }
 
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
+   //
+   // Calculate the composite object sets
+   //
+   for(mid=1; mid<=n_fcst_merge; mid++) {
 
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(first_k, -1);
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-            count++;
-         }
-      }  // end for k
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(*it, -1);
+      }
 
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+      // Increment set counter
+      collection.n_sets++;
 
-   //
-   // Done
-   //
+   } // end for mid
 
-   delete [] fcst_shape;       fcst_shape = (ShapeData *) 0;
-   delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0;
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After merging of fcst field: "
+           << fcst_split->sdebug_examine() << "\n";
+   }
 
-   mlog << Debug(4) << " After merging of fcst field: " << fcst_split->sdebug_examine() << "\n";
    return;
 }
 
+///////////////////////////////////////////////////////////////////////
+
 void ModeFuzzyEngine::do_obs_merge_thresh(const ShapeData &merge_data) {
-   int j, k, x, y;
-   int n_obs_merge, intersection;
-   int count, first_k;
-   //ShapeData obs_merge_mask;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_obs_merge;
    ShapeData obs_merge_split;
-   ShapeData * obs_shape = (ShapeData *) 0;
-   ShapeData * obs_merge_shape = (ShapeData *) 0;
 
    do_obs_splitting();
 
    //
-   // Define the forecast merge field by applying the specified threshold
-   // to the convolved field
-   //
-   // obs_merge_mask = *obs_conv;
-
-   //
-   // Threshold the forecast merge field
+   // Check dimensions
    //
-   // obs_merge_mask.threshold(conf_info.Obs->merge_thresh);
+   if(obs_split->data.nx() != merge_data.data.nx() ||
+      obs_split->data.ny() != merge_data.data.ny()) {
+      mlog << Error << "\n" << method_name
+           << "grid dimensions do not match (" << obs_split->data.nx() << ", "
+           << obs_split->data.ny() << ") != (" << merge_data.data.nx() << ", "
+           << merge_data.data.ny() << ")!\n\n";
+      exit(1);
+   }
 
    //
-   // Split up the forecast merge field
+   // Assume that the input merge data has already been split
    //
-   obs_merge_split = merge_data; //split(obs_merge_mask, n_obs_merge);
+   obs_merge_split = merge_data;
    n_obs_merge = obs_merge_split.n_objects();
 
    //
-   // Allocate space for all of the simple observation shapes and
-   // observation merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   obs_shape       = new ShapeData [n_obs];
-   obs_merge_shape = new ShapeData [n_obs_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!obs_shape || !obs_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_obs;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<obs_split->data.nxy(); j++) {
+
+      oid = obs_split->data.buf()[j];
+      mid = obs_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple observation shapes and
-   // simple merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-   }
-
-   for(j=0; j<n_obs_merge; j++) {
-      obs_merge_shape[j] = select(obs_merge_split, j+1);
+   for(oid=1; oid<=n_obs; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No observation threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object!\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_obs_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_obs; k++) {
-
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
+   for(mid=1; mid<=n_obs_merge; mid++) {
 
-         for(x=0; x<obs_merge_split.data.nx(); x++) {
-            for(y=0; y<obs_merge_split.data.ny(); y++) {
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-               if(obs_merge_shape[j].s_is_on(x, y) &&
-                  obs_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the obs object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= obs area
-         //
-         if(intersection >= obs_shape[k].area()) {
-
-            collection.make_room();
-
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
-
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(-1, first_k);
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(-1, *it);
+      }
 
-            count++;
-         }
-      }  // end for k
+      // Increment set counter
+      collection.n_sets++;
 
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+   } // end for mid
 
-   //
-   // Done
-   //
-
-   delete [] obs_shape;       obs_shape = (ShapeData *) 0;
-   delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0;
+   if(mlog.verbosity_level() >= 4) {
+      mlog << Debug(4) << " After merging of obs field: "
+           << obs_split->sdebug_examine() << "\n";
+   }
 
-   mlog << Debug(4) << " After merging of obs field: " << obs_split->sdebug_examine() << "\n";
    return;
 }
 
@@ -2114,31 +2015,22 @@ void ModeFuzzyEngine::do_obs_merge_engine(const char *default_config,
 void ModeFuzzyEngine::do_match_fcst_merge() {
    int j, k, n;
    InterestInfo junkinfo;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * obs_shape = (ShapeData *) 0;
+   ShapeData cur_shape;
 
    do_fcst_splitting();
    do_obs_splitting();
 
    clear_colors();
 
-   fcst_shape = new ShapeData [n_fcst];
-   obs_shape = new ShapeData [n_obs];
-
-   if(!fcst_shape || !obs_shape) {
-      mlog << Error << "\nModeFuzzyEngine::do_match_fcst_merge() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
-
    //
    // Do the single features
    //
    fcst_single.set_size(n_fcst);
 
    for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-      fcst_single[j].set(*fcst_raw, *fcst_thresh, fcst_shape[j],
+      cur_shape = select(*fcst_split, j+1);
+      fcst_single[j].set(*fcst_raw, *fcst_thresh,
+                         *fcst_split, cur_shape,
                          conf_info.inten_perc_value,
                          conf_info.Fcst->var_info->is_precipitation());
       fcst_single[j].object_number = j+1;
@@ -2147,8 +2039,9 @@ void ModeFuzzyEngine::do_match_fcst_merge() {
    obs_single.set_size(n_obs);
 
    for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-      obs_single[j].set(*obs_raw, *obs_thresh, obs_shape[j],
+      cur_shape = select(*obs_split, j+1);
+      obs_single[j].set(*obs_raw, *obs_thresh,
+                        *obs_split, cur_shape,
                         conf_info.inten_perc_value,
                         conf_info.Obs->var_info->is_precipitation());
       obs_single[j].object_number = j+1;
@@ -2275,9 +2168,6 @@ void ModeFuzzyEngine::do_match_fcst_merge() {
    // Done
    //
 
-   delete [] fcst_shape;  fcst_shape = (ShapeData *) 0;
-   delete [] obs_shape;   obs_shape = (ShapeData *) 0;
-
    return;
 }
 
@@ -2291,31 +2181,22 @@ void ModeFuzzyEngine::do_match_fcst_merge() {
 void ModeFuzzyEngine::do_match_only() {
    int j, k, n;
    InterestInfo junkinfo;
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * obs_shape = (ShapeData *) 0;
+   ShapeData cur_shape;
 
    do_fcst_splitting();
    do_obs_splitting();
 
    clear_colors();
 
-   fcst_shape = new ShapeData [n_fcst];
-   obs_shape = new ShapeData [n_obs];
-
-   if(!fcst_shape || !obs_shape) {
-      mlog << Error << "\nModeFuzzyEngine::do_match_only() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
-
    //
    // Do the single features
    //
    fcst_single.set_size(n_fcst);
 
    for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-      fcst_single[j].set(*fcst_raw, *fcst_thresh, fcst_shape[j],
+      cur_shape = select(*fcst_split, j+1);
+      fcst_single[j].set(*fcst_raw, *fcst_thresh,
+                         *fcst_split, cur_shape,
                          conf_info.inten_perc_value,
                          conf_info.Fcst->var_info->is_precipitation());
       fcst_single[j].object_number = j+1;
@@ -2324,8 +2205,9 @@ void ModeFuzzyEngine::do_match_only() {
    obs_single.set_size(n_obs);
 
    for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-      obs_single[j].set(*obs_raw, *obs_thresh, obs_shape[j],
+      cur_shape = select(*obs_split, j+1);
+      obs_single[j].set(*obs_raw, *obs_thresh,
+                        *obs_split, cur_shape,
                         conf_info.inten_perc_value,
                         conf_info.Obs->var_info->is_precipitation());
       obs_single[j].object_number = j+1;
@@ -2446,9 +2328,6 @@ void ModeFuzzyEngine::do_match_only() {
    // Done
    //
 
-   delete [] fcst_shape;   fcst_shape = (ShapeData *) 0;
-   delete [] obs_shape;    obs_shape = (ShapeData *) 0;
-
    return;
 }
 
@@ -2537,8 +2416,7 @@ void ModeFuzzyEngine::do_obs_clus_splitting() {
 
 void ModeFuzzyEngine::do_cluster_features() {
    int j;
-   ShapeData * fcst_clus_shape = (ShapeData *) 0;
-   ShapeData * obs_clus_shape  = (ShapeData *) 0;
+   ShapeData cur_shape;
 
    if(need_fcst_clus_split) do_fcst_clus_splitting();
    if(need_obs_clus_split)  do_obs_clus_splitting();
@@ -2548,16 +2426,6 @@ void ModeFuzzyEngine::do_cluster_features() {
    //
    n_clus = collection.n_sets;
 
-   fcst_clus_shape = new ShapeData [n_clus];
-   obs_clus_shape  = new ShapeData [n_clus];
-
-   if(!fcst_clus_shape || !obs_clus_shape) {
-
-      mlog << Error << "\nModeFuzzyEngine::do_cluster_features() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
-
    //
    // Do the single features for clusters
    //
@@ -2565,16 +2433,18 @@ void ModeFuzzyEngine::do_cluster_features() {
     obs_cluster.set_size(n_clus);
 
    for(j=0; j<n_clus; j++) {
-      fcst_clus_shape[j] = select(*fcst_clus_split, j+1);
-      fcst_cluster[j].set(*fcst_raw, *fcst_thresh, fcst_clus_shape[j],
-                       conf_info.inten_perc_value,
-                       conf_info.Fcst->var_info->is_precipitation());
+      cur_shape = select(*fcst_clus_split, j+1);
+      fcst_cluster[j].set(*fcst_raw, *fcst_thresh,
+                          *fcst_clus_split, cur_shape,
+                          conf_info.inten_perc_value,
+                          conf_info.Fcst->var_info->is_precipitation());
       fcst_cluster[j].object_number = j+1;
 
-      obs_clus_shape[j] = select(*obs_clus_split, j+1);
-      obs_cluster[j].set(*obs_raw, *obs_thresh, obs_clus_shape[j],
-                      conf_info.inten_perc_value,
-                      conf_info.Obs->var_info->is_precipitation());
+      cur_shape = select(*obs_clus_split, j+1);
+      obs_cluster[j].set(*obs_raw, *obs_thresh,
+                         *obs_clus_split, cur_shape,
+                         conf_info.inten_perc_value,
+                         conf_info.Obs->var_info->is_precipitation());
       obs_cluster[j].object_number = j+1;
    }
 
@@ -2585,7 +2455,7 @@ void ModeFuzzyEngine::do_cluster_features() {
 
    for(j=0; j<n_clus; j++) {
       pair_cluster[j].set(fcst_cluster[j], obs_cluster[j],
-                       conf_info.max_centroid_dist);
+                          conf_info.max_centroid_dist);
       pair_cluster[j].pair_number = j+1;
    }
 
@@ -2607,9 +2477,6 @@ void ModeFuzzyEngine::do_cluster_features() {
    // Done
    //
 
-   delete [] fcst_clus_shape; fcst_clus_shape = (ShapeData *) 0;
-   delete [] obs_clus_shape;  obs_clus_shape  = (ShapeData *) 0;
-
    return;
 }
 
diff --git a/src/libcode/vx_shapedata/interest.cc b/src/libcode/vx_shapedata/interest.cc
index 30cba17815..6992ed8c7f 100644
--- a/src/libcode/vx_shapedata/interest.cc
+++ b/src/libcode/vx_shapedata/interest.cc
@@ -15,9 +15,9 @@
 //
 //  Mod#   Date      Name           Description
 //  ----   ----      ----           -----------
-//  000    04-15-05  Halley Gotway
-//
-//  001    01-10-12  Bullock        Ported to new repository
+//  000    04/15/05  Halley Gotway
+//  001    01/10/12  Bullock        Ported to new repository
+//  002    11/02/23  Halley Gotway  MET #2724 improve efficiency
 //
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -128,7 +128,7 @@ void SingleFeature::clear()
 
    Raw    = (const ShapeData *) 0;
    Thresh = (const ShapeData *) 0;
-   Mask   = (const ShapeData *) 0;
+   Split  = (const ShapeData *) 0;
 
    convex_hull.clear();
 
@@ -174,7 +174,7 @@ void SingleFeature::assign(const SingleFeature & s)
 
    Raw    = s.Raw;
    Thresh = s.Thresh;
-   Mask   = s.Mask;
+   Split  = s.Split;
 
    convex_hull = s.convex_hull;
 
@@ -190,35 +190,35 @@ void SingleFeature::assign(const SingleFeature & s)
 
 ////////////////////////////////////////////////////////////////////////
 
-void SingleFeature::set(const ShapeData &raw, const ShapeData &thresh,
-                        const ShapeData &mask, const int perc,
-                        bool precip_flag)
+void SingleFeature::set(const ShapeData &raw_sd, const ShapeData &thresh_sd,
+                        const ShapeData &split_sd, const ShapeData &mask_sd,
+                        const int perc, bool precip_flag)
 
 {
 
    int i;
-   ShapeData split_wd, obj_wd;
+   ShapeData cur_split_sd, cur_obj_sd;
 
    clear();
 
-   Raw    = &raw;
-   Thresh = &thresh;
-   Mask   = &mask;
+   Raw    = &raw_sd;
+   Thresh = &thresh_sd;
+   Split  = &split_sd;
 
    //
    // Centroid
    //
-   Mask->centroid(centroid_x, centroid_y);
+   mask_sd.centroid(centroid_x, centroid_y);
 
    //
    // Axis angle
    //
-   axis_ang = Mask->angle_degrees();
+   axis_ang = mask_sd.angle_degrees();
 
    //
    // Length and Width
    //
-   Mask->calc_length_width(length, width);
+   mask_sd.calc_length_width(length, width);
 
    //
    // Aspect ratio
@@ -228,28 +228,28 @@ void SingleFeature::set(const ShapeData &raw, const ShapeData &thresh,
    //
    // Object area
    //
-   area = Mask->area();
+   area = mask_sd.area();
 
    //
    // Object threshold area: the area of the raw field inside the mask
    // area that meets the threshold criteria
    //
-   area_thresh = (double) ShapeData_intersection(*Thresh, *Mask);
+   area_thresh = (double) ShapeData_intersection(*Thresh, mask_sd);
 
    //
    // Curvature, Curvature_x, Curvature_y
    //
-   curvature = Mask->curvature(curvature_x, curvature_y);
+   curvature = mask_sd.curvature(curvature_x, curvature_y);
 
    //
    // Complexity
    //
-   complexity = Mask->complexity();
+   complexity = mask_sd.complexity();
 
    //
    // Compute the Intensity Percentiles
    //
-   get_percentiles(intensity_ptile, raw, mask, perc, precip_flag);
+   get_percentiles(intensity_ptile, raw_sd, mask_sd, perc, precip_flag);
 
    //
    // User Percentile
@@ -259,17 +259,17 @@ void SingleFeature::set(const ShapeData &raw, const ShapeData &thresh,
    //
    // Convex hull
    //
-   convex_hull = Mask->convex_hull();
+   convex_hull = mask_sd.convex_hull();
 
    //
    // Boundary:
    // Split the mask field and store the boundary for each object.
    //
-   split_wd = split(mask, n_bdy);
+   cur_split_sd = split(mask_sd, n_bdy);
    boundary = new Polyline [n_bdy];
    for(i=0; i<n_bdy; i++) {
-      obj_wd      = select(split_wd, i+1);
-      boundary[i] = obj_wd.single_boundary();
+      cur_obj_sd  = select(cur_split_sd, i+1);
+      boundary[i] = cur_obj_sd.single_boundary();
    }
 
    //
@@ -407,7 +407,6 @@ void PairFeature::set(const SingleFeature &fcst,
    Fcst = &fcst;
    Obs  = &obs;
 
-   int x, y;
    int fcst_on, obs_on;
    double dx, dy;
    double a1, a2;
@@ -461,23 +460,21 @@ void PairFeature::set(const SingleFeature &fcst,
    area_ratio    = min( (Obs->area)/(Fcst->area),
                         (Fcst->area)/(Obs->area) );
 
-
    //
    // Intersection, union, and symmetric diff areas
    //
    intersection_area = union_area = 0.0;
    symmetric_diff = 0.0;
-   for(x=0; x<(Fcst->Mask->data.nx()); ++x) {
-      for(y=0; y<(Fcst->Mask->data.ny()); ++y) {
+   int nxy = Fcst->Split->data.nxy();
+   for(i=0; i<nxy; ++i) {
 
-         fcst_on = Fcst->Mask->s_is_on(x, y);
-         obs_on  =  Obs->Mask->s_is_on(x, y);
+      fcst_on = (nint(Fcst->Split->data.data()[i]) == Fcst->object_number ? 1 : 0);
+      obs_on  = (nint( Obs->Split->data.data()[i]) ==  Obs->object_number ? 1 : 0);
 
-         if(fcst_on && obs_on) intersection_area++;
-         if(fcst_on || obs_on) union_area++;
-         if((fcst_on && !obs_on) ||
-            (!fcst_on && obs_on)) symmetric_diff++;
-      }
+      if(   fcst_on &&  obs_on)   intersection_area++;
+      if(   fcst_on ||  obs_on)   union_area++;
+      if( ( fcst_on && !obs_on) ||
+          (!fcst_on &&  obs_on) ) symmetric_diff++;
    }
 
    //
diff --git a/src/libcode/vx_shapedata/interest.h b/src/libcode/vx_shapedata/interest.h
index de8832ed2b..6fd6f9cc2d 100644
--- a/src/libcode/vx_shapedata/interest.h
+++ b/src/libcode/vx_shapedata/interest.h
@@ -49,13 +49,13 @@ class SingleFeature {
 
       void clear();
 
-      void set(const ShapeData &raw,  const ShapeData &thresh,
-               const ShapeData &mask, const int perc,
-               const bool precip_flag);
+      void set(const ShapeData &raw, const ShapeData &thresh,
+               const ShapeData &split, const ShapeData &mask,
+               const int perc, const bool precip_flag);
 
       const ShapeData * Raw;    //  NOT allocated, so DON'T delete!
       const ShapeData * Thresh; //  NOT allocated, so DON'T delete!
-      const ShapeData * Mask;   //  NOT allocated, so DON'T delete!
+      const ShapeData * Split;  //  NOT allocated, so DON'T delete!
 
       int object_number;
 
diff --git a/src/libcode/vx_shapedata/shapedata.cc b/src/libcode/vx_shapedata/shapedata.cc
index 9051d00819..88b63e195b 100644
--- a/src/libcode/vx_shapedata/shapedata.cc
+++ b/src/libcode/vx_shapedata/shapedata.cc
@@ -15,8 +15,9 @@
 //
 //   Mod#   Date      Name           Description
 //   ----   ----      ----           -----------
-//   000    11-05-31  Halley Gotway  Adapated from wrfdata.cc.
-//   001    14-05-29  Halley Gotway  Add ShapeData::n_objects()
+//   000    05/31/11  Halley Gotway  Adapated from wrfdata.cc
+//   001    05/29/14  Halley Gotway  Add ShapeData::n_objects()
+//   002    11/02/23  Halley Gotway  MET #2724 add OpenMP to convolution
 //
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -513,160 +514,117 @@ double ShapeData::get_attr(const ConcatString &attr_name,
 
 ///////////////////////////////////////////////////////////////////////////////
 
+void ShapeData::conv_filter_circ(int diameter, double vld_thresh) {
+   const char *method_name = "ShapeData::conv_filter_circ() -> ";
+   GridPoint *gp = nullptr;
+   int x, y, count, n_vld;
+   double v, v_sum;
+   DataPlane conv_dp;
 
-void ShapeData::conv_filter_circ(int diameter, double vld_thresh)
-
-{
-
-int x, y, xx, yy, u, v;
-int dn, fn, nn;
-int vpr, upr;
-int count, bd_count;
-double center, cur, sum;
-double dx, dy, dist;
-double vld_ratio;
-const int nx = data.nx();
-const int ny = data.ny();
-bool * f = (bool *) 0;
-bool center_bad = false;
-DataPlane in_data = data;
-const bool vld_thresh_one = is_eq(vld_thresh, 1.0);
-
-
-if ( (diameter%2 == 0) || (diameter < 3) )  {
-
-   mlog << Error << "\nShapeData::conv_filter_circ() -> "
-        << "diameter must be odd and >= 3 ... diameter = "
-        << diameter << "\n\n";
-
-   exit(1);
-
-}
-
-const int radius = (diameter - 1)/2;
-
-const vector<double> * in  = &(in_data.Data);
-      vector<double> * out = &(data.Data);
-
-f = new bool [diameter*diameter];
-
-   //
-   //  set up the filter
-   //
-
-for (y=0; y<diameter; ++y)  {
-
-   dy = y - radius;
-
-   for (x=0; x<diameter; ++x)  {
-
-      dx = x - radius;
-
-      dist = sqrt( dx*dx + dy*dy );
-
-      fn = STANDARD_XY_YO_N(diameter, x, y) ;
-
-      f[fn] = (dist <= radius);
-
+   // Check the diameter
+   if(diameter%2 == 0 || diameter < 3) {
+      mlog << Error << "\n" << method_name
+           << "diameter must be odd and >= 3 ... diameter = "
+           << diameter << "\n\n";
+      exit(1);
    }
 
-}
-
-   //
-   //  do the convolution
-   //
-
-dn = -1;
-
-for(y=0; y<ny; y++) {
-
-   for(x=0; x<nx; x++) {
-
-      ++dn;
-
-         //
-         // If the bad data threshold is set to zero and the center of the
-         // convolution radius contains bad data, set the convolved value to
-         // bad data and continue.
-         //
-
-      center = (*in)[dn];
-
-      center_bad = ::is_bad_data(center);
-
-      if ( center_bad && vld_thresh_one ) { (*out)[dn] = bad_data_double;  continue; }
-
-      sum      = 0.0;
-      count    = 0;
-      bd_count = 0;
-
-      for (v=-radius; v<=radius; ++v) {
-
-         yy = y + v;
-
-         if ( (yy < 0) || (yy >= ny) )  continue;
-
-         vpr = v + radius;
-
-         for(u=-radius; u<=radius; ++u) {
-
-            xx = x + u;
-
-            if ( (xx < 0) || (xx >= nx) )  continue;
-
-            upr = u + radius;
-
-            fn = STANDARD_XY_YO_N(diameter, upr, vpr);
-
-            if ( !(f[fn]) )  continue;
-
-            nn = STANDARD_XY_YO_N(nx, xx, yy) ;
-
-            cur = (*in)[nn];
-
-            if( ::is_bad_data(cur) ) { bd_count++;  continue; }
-
-            sum += cur;
-
-            count++;
-
-         } // for v
-
-      } // for u
-
-         //
-         //  If the center of the convolution contains bad data and the ratio
-         //  of bad data in the convolution area is too high, set the convoled
-         //  value to bad data.
-         //
-
-      if ( count == 0 )  sum = bad_data_double;
-      else {
+#pragma omp parallel default(none)                   \
+   shared(mlog, data, conv_dp, diameter, vld_thresh) \
+   private(x, y, count, n_vld, v, v_sum, gp)
+   {
+
+      // Build the grid template with shape circle and wrap_lon false
+      GridTemplateFactory gtf;
+      GridTemplate* gt = gtf.buildGT(GridTemplateFactory::GridTemplate_Circle,
+                                     diameter, false);
+
+#pragma omp single
+      {
+         // Initialize the convolved field to bad data
+         conv_dp = data;
+         conv_dp.set_constant(bad_data_double);
+      }
 
-         vld_ratio = ((double) count)/(bd_count + count);
+      // Compute the convolved values
+#pragma omp for schedule (static)
+      for(x=0; x<data.nx(); x++) {
+         for(y=0; y<data.ny(); y++) {
+
+            // For a new column, reset the grid template and counts
+            if(y == 0) {
+
+               // Initialize counts and sum
+               count = 0;
+               n_vld = 0;
+               v_sum = 0.0;
+
+               // Sum all the points
+               for(gp  = gt->getFirstInGrid(x, y, data.nx(), data.ny());
+                   gp != nullptr;
+                   gp  = gt->getNextInGrid()) {
+                  v = data.get(gp->x, gp->y);
+                  count += 1;
+                  if(::is_bad_data(v)) continue;
+                  n_vld += 1;
+                  v_sum += v;
+               }
+            }
+            // Subtract off the bottom edge, shift up, and add the top
+            else {
+
+               // Subtract points from the the bottom edge
+               for(gp  = gt->getFirstInBotEdge();
+                   gp != nullptr;
+                   gp  = gt->getNextInBotEdge()) {
+                  v = data.get(gp->x, gp->y);
+                  count -= 1;
+                  if(::is_bad_data(v)) continue;
+                  n_vld -= 1;
+                  v_sum -= v;
+               }
+
+               // Increment Y
+               gt->incBaseY(1);
+
+               // Add points from the the top edge
+               for(gp  = gt->getFirstInTopEdge();
+                   gp != nullptr;
+                   gp  = gt->getNextInTopEdge()) {
+                  v = data.get(gp->x, gp->y);
+                  count += 1;
+                  if(::is_bad_data(v)) continue;
+                  n_vld += 1;
+                  v_sum += v;
+               }
+            }
 
-         if ( center_bad && (vld_ratio < vld_thresh) )  sum = bad_data_double;
-         else                                           sum /= count;
+            //  If the center of the convolution contains bad data and the ratio
+            //  of bad data in the convolution area is too high, set the convoled
+            //  value to bad data.
+            if(count == 0 || n_vld == 0)                v = bad_data_double;
+            else if(::is_bad_data(data.get(x, y)) &&
+                    (double)(n_vld)/count < vld_thresh) v = bad_data_double;
+            else                                        v = (double) v_sum/n_vld;
+            conv_dp.set(v, x, y);
 
-      }
+         } // end for y
 
-      (*out)[dn] = sum;
+         // Increment X
+         if(x < (data.nx() - 1)) gt->incBaseX(1);
 
-   } // for y
+      } // end for x
 
-} // for x
+      delete gt;
 
-   //
-   //  done
-   //
+   } // End of omp parallel
 
-if ( f )  { delete [] f;   f = (bool *) 0; }
-
-return;
+   // Save the result
+   data = conv_dp;
 
+   return;
 }
 
-
 ////////////////////////////////////////////////////////////////////////
 
 
diff --git a/src/tools/core/ensemble_stat/ensemble_stat.cc b/src/tools/core/ensemble_stat/ensemble_stat.cc
index 62885a5f5e..b778e97f2f 100644
--- a/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -99,8 +99,6 @@ using namespace netCDF;
 #include "nc_obs_util.h"
 #include "nc_point_obs_in.h"
 
-#include "handle_openmp.h"
-
 #ifdef WITH_PYTHON
 #include "data2d_nc_met.h"
 #include "pointdata_python.h"
@@ -177,9 +175,6 @@ static void set_compress(const StringArray &);
 
 int met_main(int argc, char *argv[]) {
 
-   // Set up OpenMP (if enabled)
-   init_openmp();
-
    // Process the command line arguments
    process_command_line(argc, argv);
 
diff --git a/src/tools/core/mode/mode.cc b/src/tools/core/mode/mode.cc
index 88e18d0b40..a9277a0884 100644
--- a/src/tools/core/mode/mode.cc
+++ b/src/tools/core/mode/mode.cc
@@ -56,6 +56,7 @@
 //   019    04/01/19  Fillmore       Add FCST and OBS units.
 //   020    07/06/22  Howard Soh     METplus-Internal #19 Rename main to met_main
 //   021    06/09/23  Albo           Major changes for multivariate mode
+//   022    11/02/23  Halley Gotway  MET #2724 add OpenMP to convolution
 //
 ////////////////////////////////////////////////////////////////////////
 
@@ -73,6 +74,7 @@ using namespace std;
 #include <sys/types.h>
 
 #include "main.h"
+#include "handle_openmp.h"
 #include "string_array.h"
 #include "mode_usage.h"
 #include "mode_frontend.h"
@@ -131,6 +133,9 @@ int met_main(int argc, char * argv [])
    string s;
    const char * user_config_filename = 0;
 
+   // Set up OpenMP (if enabled)
+   init_openmp();
+
    for (j=0,n=0; j<argc; ++j)  {
 
       //