From 3b0d6d27153aa9230e787e5c72cbc297219b812b Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 8 Oct 2021 11:38:43 -0600
Subject: [PATCH 01/19] Per #1905, no real change. Just spacing

---
 met/src/tools/other/gen_ens_prod/gen_ens_prod.cc | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc b/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc
index ce956467af..4fc7491c86 100644
--- a/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc
+++ b/met/src/tools/other/gen_ens_prod/gen_ens_prod.cc
@@ -904,8 +904,7 @@ void write_ens_var_int(int i_var, int *ens_data, const DataPlane &dp,
 ////////////////////////////////////////////////////////////////////////
 
 void write_ens_data_plane(int i_var, const DataPlane &ens_dp, const DataPlane &dp,
-                         const char *type_str,
-                         const char *long_name_str) {
+                          const char *type_str, const char *long_name_str) {
 
    // Allocate memory for this data
    float *ens_data = new float [nxy];

From 144ac67fcaabeb5d7ccd85df6ecd6c13cc1f50b1 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 8 Oct 2021 11:44:55 -0600
Subject: [PATCH 02/19] Per #1905, add the -ctrl command line option to
 ensemble_stat and update the logic to use it the same way gen_ens_prod does.

---
 .../tools/core/ensemble_stat/ensemble_stat.cc | 443 +++++++++---------
 .../tools/core/ensemble_stat/ensemble_stat.h  |  17 +-
 2 files changed, 232 insertions(+), 228 deletions(-)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
index 6dfc98dee3..77420adca5 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -61,6 +61,7 @@
 //                    statistics in the ECNT line type.
 //   031    09/13/21  Seth Linden    Changed obs_qty to obs_qty_inc.
 //                    Added code for obs_qty_exc.
+//   032    10/07/21  Halley Gotway  MET #1905 Add -ctrl option.
 //
 ////////////////////////////////////////////////////////////////////////
 
@@ -124,7 +125,7 @@ static void do_rps               (const EnsembleStatVxOpt &,
                                   const PairDataEnsemble *);
 
 static void clear_counts();
-static void track_counts(int, const DataPlane &);
+static void track_counts(int, const DataPlane &, bool);
 
 static ConcatString get_ens_mn_var_name(int);
 
@@ -135,10 +136,12 @@ static void setup_table     (AsciiTable &);
 static void build_outfile_name(unixtime, const char *,
                                ConcatString &);
 static void write_ens_nc(int, DataPlane &);
-static void write_ens_var_float(int, float *, DataPlane &,
+static void write_ens_var_float(int, float *, const DataPlane &,
                                 const char *, const char *);
-static void write_ens_var_int(int, int *, DataPlane &,
+static void write_ens_var_int(int, int *, const DataPlane &,
                               const char *, const char *);
+static void write_ens_data_plane(int, const DataPlane &, const DataPlane &,
+                                 const char *, const char *);
 
 static void write_orank_nc(PairDataEnsemble &, DataPlane &, int, int, int);
 static void write_orank_var_float(int, int, int, float *, DataPlane &,
@@ -146,8 +149,8 @@ static void write_orank_var_float(int, int, int, float *, DataPlane &,
 static void write_orank_var_int(int, int, int, int *, DataPlane &,
                                 const char *, const char *);
 
-static void add_var_att_local(VarInfo *, NcVar *, bool is_int, DataPlane &,
-                              const char *, const char *);
+static void add_var_att_local(VarInfo *, NcVar *, bool is_int,
+                              const DataPlane &, const char *, const char *);
 
 static void finish_txt_files();
 static void clean_up();
@@ -157,6 +160,7 @@ static void usage();
 static void set_grid_obs(const StringArray &);
 static void set_point_obs(const StringArray &);
 static void set_ens_mean(const StringArray & a);
+static void set_ctrl_file(const StringArray &);
 static void set_obs_valid_beg(const StringArray &);
 static void set_obs_valid_end(const StringArray &);
 static void set_outdir(const StringArray &);
@@ -220,6 +224,7 @@ void process_command_line(int argc, char **argv) {
    cline.add(set_grid_obs,      "-grid_obs", 1);
    cline.add(set_point_obs,     "-point_obs", 1);
    cline.add(set_ens_mean,      "-ens_mean", 1);
+   cline.add(set_ctrl_file,     "-ctrl", 1);
    cline.add(set_obs_valid_beg, "-obs_valid_beg", 1);
    cline.add(set_obs_valid_end, "-obs_valid_end", 1);
    cline.add(set_outdir,        "-outdir", 1);
@@ -370,11 +375,11 @@ void process_command_line(int argc, char **argv) {
    shc.set_model(conf_info.model.c_str());
 
    // Allocate arrays to store threshold counts
-   thresh_count_na       = new NumArray   [conf_info.get_max_n_thresh()];
-   thresh_nbrhd_count_na = new NumArray * [conf_info.get_max_n_thresh()];
+   thresh_cnt_na       = new NumArray   [conf_info.get_max_n_thresh()];
+   thresh_nbrhd_cnt_na = new NumArray * [conf_info.get_max_n_thresh()];
 
    for(i=0; i<conf_info.get_max_n_thresh(); i++) {
-      thresh_nbrhd_count_na[i] = new NumArray [conf_info.get_n_nbrhd()];
+      thresh_nbrhd_cnt_na[i] = new NumArray [conf_info.get_n_nbrhd()];
    }
 
    // List the input ensemble files
@@ -384,6 +389,10 @@ void process_command_line(int argc, char **argv) {
       mlog << "   " << ens_file_list[i]  << "\n";
    }
 
+   // List the control member file
+   mlog << Debug(1) << "Ensemble Control: "
+        << (ctrl_file.empty() ? "None" : ctrl_file.c_str()) << "\n";
+
    // List the input gridded observations files
    if(grid_obs_file_list.n() > 0) {
       mlog << Debug(1) << "Gridded Observation Files["
@@ -755,8 +764,8 @@ bool get_data_plane_array(const char *infile, GrdFileType ftype,
 
 void process_ensemble() {
    int i, j;
-   bool reset;
-   DataPlane ens_dp;
+   bool need_reset;
+   DataPlane ens_dp, ctrl_dp;
    unixtime max_init_ut = bad_data_ll;
 
    // Loop through each of the ensemble fields to be processed
@@ -766,8 +775,11 @@ void process_ensemble() {
            << "Processing ensemble field: "
            << conf_info.ens_info[i]->magic_str() << "\n";
 
+      // Need to reinitialize counts and sums for each ensemble field
+      need_reset = true;
+
       // Loop through each of the input forecast files
-      for(j=0, reset=true; j<ens_file_list.n(); j++) {
+      for(j=0; j<ens_file_list.n(); j++) {
 
          // Skip bad data files
          if(!ens_file_vld[j]) continue;
@@ -782,13 +794,40 @@ void process_ensemble() {
          }
 
          // Reset the running sums and counts
-         if(reset) {
+         if(need_reset) {
+
+            need_reset = false;
+
+            // Reset the running sums and counts
             clear_counts();
-            reset = false;
-         }
+
+            // Read ensemble control member data, if provided
+            if(ctrl_file.nonempty()) {
+
+               // Error out if missing
+               if(!get_data_plane(ctrl_file.c_str(), etype,
+                                  conf_info.ens_info[i], ctrl_dp, true)) {
+                  mlog << Error << "\nprocess_ensemble() -> "
+                       << "control member ensemble field \""
+                       << conf_info.ens_info[i]->magic_str()
+                       << "\" not found in file \"" << ctrl_file << "\"\n\n";
+                  exit(1);
+               }
+
+               // Apply current data to the running sums and counts
+               track_counts(i, ctrl_dp, true);
+
+            } // end if ctrl_file
+
+            mlog << Debug(3)
+                 << "Found " << (ctrl_dp.is_empty() ? 0 : 1)
+                 << " control member field(s) for \""
+                 << conf_info.ens_info[i]->magic_str() << "\".\n";
+
+         } // end if need_reset
 
          // Apply current data to the running sums and counts
-         track_counts(i, ens_dp);
+         track_counts(i, ens_dp, false);
 
          // Keep track of the maximum initialization time
          if(is_bad_data(max_init_ut) || ens_dp.init() > max_init_ut) {
@@ -1813,7 +1852,7 @@ void process_grid_scores(int i_vx,
    ObsErrorEntry *e = (ObsErrorEntry *) 0;
 
    // Allocate memory in one big chunk based on grid size
-   pd.extend(grid.nx()*grid.ny());
+   pd.extend(nxy);
 
    // Climatology flags
    bool emn_flag = (emn_dp.nx() == obs_dp.nx() &&
@@ -1955,46 +1994,21 @@ void do_rps(const EnsembleStatVxOpt &vx_opt,
 ////////////////////////////////////////////////////////////////////////
 
 void clear_counts() {
-   int i, j, k;
+   int i, j;
 
-   // Allocate memory in one big chunk based on grid size, if needed
-   count_na.extend(nxy);
-   min_na.extend(nxy);
-   max_na.extend(nxy);
-   sum_na.extend(nxy);
-   sum_sq_na.extend(nxy);
-   for(i=0; i<conf_info.get_max_n_thresh(); i++) {
-      thresh_count_na[i].extend(nxy);
-      for(j=0; j<conf_info.get_n_nbrhd(); j++) {
-         thresh_nbrhd_count_na[i][j].extend(nxy);
-      }
-   }
+   cnt_na.set_const(0.0, nxy);
+   min_na.set_const(bad_data_double, nxy);
+   max_na.set_const(bad_data_double, nxy);
+   sum_na.set_const(0.0, nxy);
+
+   stdev_cnt_na.set_const(0.0, nxy);
+   stdev_sum_na.set_const(0.0, nxy);
+   stdev_ssq_na.set_const(0.0, nxy);
 
-   // Erase existing values
-   count_na.erase();
-   min_na.erase();
-   max_na.erase();
-   sum_na.erase();
-   sum_sq_na.erase();
    for(i=0; i<conf_info.get_max_n_thresh(); i++) {
-      thresh_count_na[i].erase();
+      thresh_cnt_na[i].set_const(0.0, nxy);
       for(j=0; j<conf_info.get_n_nbrhd(); j++) {
-         thresh_nbrhd_count_na[i][j].erase();
-      }
-   }
-
-   // Initialize arrays
-   for(i=0; i<nxy; i++) {
-      count_na.add(0);
-      min_na.add(bad_data_double);
-      max_na.add(bad_data_double);
-      sum_na.add(0.0);
-      sum_sq_na.add(0.0);
-      for(j=0; j<conf_info.get_max_n_thresh(); j++) {
-         thresh_count_na[j].add(0);
-         for(k=0; k<conf_info.get_n_nbrhd(); k++) {
-            thresh_nbrhd_count_na[j][k].add(0);
-         }
+         thresh_nbrhd_cnt_na[i][j].set_const(0.0, nxy);
       }
    }
 
@@ -2003,47 +2017,48 @@ void clear_counts() {
 
 ////////////////////////////////////////////////////////////////////////
 
-void track_counts(int i_vx, const DataPlane &dp) {
+void track_counts(int i_vx, const DataPlane &ens_dp, bool is_ctrl) {
    int i, j, k;
-   double v;
-
-   // Pointers to data buffers for faster access
-   const double *Data = dp.data();
-   double *CountBuf   = count_na.buf();
-   double *MinBuf     = min_na.buf();
-   double *MaxBuf     = max_na.buf();
-   double *SumBuf     = sum_na.buf();
-   double *SumSqBuf   = sum_sq_na.buf();
+   double ens;
 
    // Ensemble thresholds
-   const int Nthresh = conf_info.ens_ta[i_vx].n();
-   SingleThresh *ThreshBuf = conf_info.ens_ta[i_vx].buf();
+   const int n_thr = conf_info.ens_ta[i_vx].n();
+   SingleThresh *thr_buf = conf_info.ens_ta[i_vx].buf();
 
    // Increment counts for each grid point
    for(i=0; i<nxy; i++) {
 
-      // Get current value
-      v = *Data++;
+      // Get current values
+      ens = ens_dp.data()[i];
 
-      // Skip the bad data value
-      if(is_bad_data(v)) continue;
+      // Skip bad data values
+      if(is_bad_data(ens)) continue;
 
       // Otherwise, update counts and sums
       else {
 
-         CountBuf[i] += 1;
+         // Valid data count
+         cnt_na.buf()[i] += 1;
 
-         if(v <= MinBuf[i] || is_bad_data(MinBuf[i])) MinBuf[i] = v;
-         if(v >= MaxBuf[i] || is_bad_data(MaxBuf[i])) MaxBuf[i] = v;
+         // Ensemble sum
+         sum_na.buf()[i] += ens;
 
-         SumBuf[i]   += v;
-         SumSqBuf[i] += v*v;
+         // Ensemble min and max
+         if(ens <= min_na.buf()[i] || is_bad_data(min_na.buf()[i])) min_na.buf()[i] = ens;
+         if(ens >= max_na.buf()[i] || is_bad_data(max_na.buf()[i])) max_na.buf()[i] = ens;
 
-         for(j=0; j<Nthresh; j++) {
-            if(ThreshBuf[j].check(v)) thresh_count_na[j].inc(i, 1);
+         // Standard deviation sum, sum of squares, and count, excluding control member
+         if(!is_ctrl) {
+            stdev_sum_na.buf()[i] += ens;
+            stdev_ssq_na.buf()[i] += ens*ens;
+            stdev_cnt_na.buf()[i] += 1;
          }
-      } // end else
 
+         // Event frequency
+         for(j=0; j<n_thr; j++) {
+            if(thr_buf[j].check(ens)) thresh_cnt_na[j].inc(i, 1);
+         }
+      } // end else
    } // end for i
 
    // Increment NMEP count anywhere fractional coverage > 0
@@ -2051,23 +2066,21 @@ void track_counts(int i_vx, const DataPlane &dp) {
       DataPlane frac_dp;
 
       // Loop over thresholds
-      for(i=0; i<Nthresh; i++) {
+      for(i=0; i<n_thr; i++) {
 
          // Loop over neighborhood sizes
          for(j=0; j<conf_info.get_n_nbrhd(); j++) {
 
             // Compute fractional coverage
-            fractional_coverage(dp, frac_dp,
+            fractional_coverage(ens_dp, frac_dp,
                conf_info.nbrhd_prob.width[j],
-               conf_info.nbrhd_prob.shape,
-               grid.wrap_lon(), ThreshBuf[i],
-               (const DataPlane *) 0, (const DataPlane *) 0,
+               conf_info.nbrhd_prob.shape, grid.wrap_lon(),
+               thr_buf[i], 0, 0,
                conf_info.nbrhd_prob.vld_thresh);
 
             // Increment counts
-            const double *Frac = frac_dp.data();
             for(k=0; k<nxy; k++) {
-               if(Frac[k] > 0) thresh_nbrhd_count_na[i][j].inc(k, 1);
+               if(frac_dp.data()[k] > 0) thresh_nbrhd_cnt_na[i][j].inc(k, 1);
             } // end for k
 
          } // end for j
@@ -2332,45 +2345,33 @@ void build_outfile_name(unixtime ut, const char *suffix, ConcatString &str) {
 
 ////////////////////////////////////////////////////////////////////////
 
-void write_ens_nc(int i_ens, DataPlane &dp) {
+void write_ens_nc(int i_var, DataPlane &ens_dp) {
    int i, j, k, l;
    double t, v;
    char type_str[max_str_len];
-   DataPlane prob_dp, smooth_dp;
-
-   // Arrays for storing ensemble data
-   float *ens_mean  = (float *) 0;
-   float *ens_stdev = (float *) 0;
-   float *ens_minus = (float *) 0;
-   float *ens_plus  = (float *) 0;
-   float *ens_min   = (float *) 0;
-   float *ens_max   = (float *) 0;
-   float *ens_range = (float *) 0;
-   int   *ens_vld   = (int   *) 0;
-   float *ens_prob  = (float *) 0;
+   DataPlane prob_dp, nbrhd_dp;
 
    // Allocate memory for storing ensemble data
-   ens_mean  = new float [grid.nx()*grid.ny()];
-   ens_stdev = new float [grid.nx()*grid.ny()];
-   ens_minus = new float [grid.nx()*grid.ny()];
-   ens_plus  = new float [grid.nx()*grid.ny()];
-   ens_min   = new float [grid.nx()*grid.ny()];
-   ens_max   = new float [grid.nx()*grid.ny()];
-   ens_range = new float [grid.nx()*grid.ny()];
-   ens_vld   = new int   [grid.nx()*grid.ny()];
-   ens_prob  = new float [grid.nx()*grid.ny()];
+   float * ens_mean  = new float [nxy];
+   float * ens_stdev = new float [nxy];
+   float * ens_minus = new float [nxy];
+   float * ens_plus  = new float [nxy];
+   float * ens_min   = new float [nxy];
+   float * ens_max   = new float [nxy];
+   float * ens_range = new float [nxy];
+   int   * ens_vld   = new int   [nxy];
 
    // Store the threshold for the ratio of valid data points
    t = conf_info.vld_data_thresh;
 
    // Store the data
-   for(i=0; i<count_na.n(); i++) {
+   for(i=0; i<cnt_na.n(); i++) {
 
       // Valid data count
-      ens_vld[i] = nint(count_na[i]);
+      ens_vld[i] = nint(cnt_na[i]);
 
       // Check for too much missing data
-      if((double) (count_na[i]/n_ens_vld[i_ens]) < t) {
+      if((double) (cnt_na[i]/n_ens_vld[i_var]) < t) {
          ens_mean[i]  = bad_data_float;
          ens_stdev[i] = bad_data_float;
          ens_minus[i] = bad_data_float;
@@ -2382,139 +2383,121 @@ void write_ens_nc(int i_ens, DataPlane &dp) {
       else {
 
          // Compute ensemble summary
-         ens_mean[i]  = (float) (sum_na[i]/count_na[i]);
-         ens_stdev[i] = (float) compute_stdev(sum_na[i], sum_sq_na[i], ens_vld[i]);
+         ens_mean[i]  = (float) (sum_na[i]/cnt_na[i]);
+         ens_stdev[i] = (float) compute_stdev(stdev_sum_na[i], stdev_ssq_na[i],
+                                              nint(stdev_cnt_na[i]));
          ens_minus[i] = (float) ens_mean[i] - ens_stdev[i];
          ens_plus[i]  = (float) ens_mean[i] + ens_stdev[i];
          ens_min[i]   = (float) min_na[i];
          ens_max[i]   = (float) max_na[i];
          v = max_na[i] - min_na[i];
-         if(is_eq(v, 0.0)) v = 0;
-         ens_range[i] = (float) v;
+         ens_range[i] = (is_eq(max_na[i], min_na[i]) ?
+                         0.0 : (float) max_na[i] - min_na[i]);
       }
    } // end for i
 
    // Add the ensemble mean, if requested
    if(conf_info.nc_info.do_mean) {
-      write_ens_var_float(i_ens, ens_mean, dp,
+      write_ens_var_float(i_var, ens_mean, ens_dp,
                           "ENS_MEAN",
                           "Ensemble Mean");
    }
 
    // Add the ensemble standard deviation, if requested
    if(conf_info.nc_info.do_stdev) {
-      write_ens_var_float(i_ens, ens_stdev, dp,
+      write_ens_var_float(i_var, ens_stdev, ens_dp,
                           "ENS_STDEV",
                           "Ensemble Standard Deviation");
    }
 
    // Add the ensemble mean minus one standard deviation, if requested
    if(conf_info.nc_info.do_minus) {
-      write_ens_var_float(i_ens, ens_minus, dp,
+      write_ens_var_float(i_var, ens_minus, ens_dp,
                           "ENS_MINUS",
                           "Ensemble Mean Minus 1 Standard Deviation");
    }
 
    // Add the ensemble mean plus one standard deviation, if requested
    if(conf_info.nc_info.do_plus) {
-      write_ens_var_float(i_ens, ens_plus, dp,
+      write_ens_var_float(i_var, ens_plus, ens_dp,
                           "ENS_PLUS",
                           "Ensemble Mean Plus 1 Standard Deviation");
    }
 
    // Add the ensemble minimum value, if requested
    if(conf_info.nc_info.do_min) {
-      write_ens_var_float(i_ens, ens_min, dp,
+      write_ens_var_float(i_var, ens_min, ens_dp,
                           "ENS_MIN",
                           "Ensemble Minimum");
    }
 
    // Add the ensemble maximum value, if requested
    if(conf_info.nc_info.do_max) {
-      write_ens_var_float(i_ens, ens_max, dp,
+      write_ens_var_float(i_var, ens_max, ens_dp,
                           "ENS_MAX",
                           "Ensemble Maximum");
    }
 
    // Add the ensemble range, if requested
    if(conf_info.nc_info.do_range) {
-      write_ens_var_float(i_ens, ens_range, dp,
+      write_ens_var_float(i_var, ens_range, ens_dp,
                           "ENS_RANGE",
                           "Ensemble Range");
    }
 
    // Add the ensemble valid data count, if requested
    if(conf_info.nc_info.do_vld) {
-      write_ens_var_int(i_ens, ens_vld, dp,
+      write_ens_var_int(i_var, ens_vld, ens_dp,
                         "ENS_VLD",
                         "Ensemble Valid Data Count");
    }
 
    // Add the ensemble relative frequencies and neighborhood probabilities, if requested
-   if(conf_info.nc_info.do_freq || conf_info.nc_info.do_nep) {
+   if(conf_info.nc_info.do_freq ||
+      conf_info.nc_info.do_nep) {
 
-      prob_dp.set_size(dp.nx(), dp.ny());
+      prob_dp.set_size(grid.nx(), grid.ny());
 
       // Loop through each threshold
-      for(i=0; i<conf_info.ens_ta[i_ens].n(); i++) {
+      for(i=0; i<conf_info.ens_ta[i_var].n(); i++) {
 
          // Initialize
          prob_dp.erase();
 
-         // Output variable name
-         snprintf(type_str, sizeof(type_str), "ENS_FREQ_%s",
-                  conf_info.ens_ta[i_ens][i].get_abbr_str().contents().c_str());
-
          // Compute the ensemble relative frequency
-         for(j=0; j<count_na.n(); j++) {
-
-            // Check for too much missing data
-            if((double) (count_na[j]/n_ens_vld[i_ens]) < t) {
-               ens_prob[j] = bad_data_float;
-            }
-            else {
-               ens_prob[j] = (float) (thresh_count_na[i][j]/count_na[j]);
-            }
-
-            // Also store value in a DataPlane object
-            prob_dp.buf()[j] = ens_prob[j];
-
-         } // end for j
+         for(j=0; j<cnt_na.n(); j++) {
+            prob_dp.buf()[j] = ((double) (cnt_na[j]/n_ens_vld[i_var]) < t ?
+                                bad_data_double :
+                                (double) (thresh_cnt_na[i][j]/cnt_na[j]));
+         }
 
-         // Write the ensemble relative frequency
+         // Write ensemble relative frequency
          if(conf_info.nc_info.do_freq) {
-            write_ens_var_float(i_ens, ens_prob, dp,
-                                type_str,
-                                "Ensemble Relative Frequency");
+            snprintf(type_str, sizeof(type_str), "ENS_FREQ_%s",
+                     conf_info.ens_ta[i_var][i].get_abbr_str().contents().c_str());
+            write_ens_data_plane(i_var, prob_dp, ens_dp, type_str,
+                                 "Ensemble Relative Frequency");
          }
 
          // Write the neighborhood ensemble probability
          if(conf_info.nc_info.do_nep) {
-
             GaussianInfo info;
-            // Compute the mean neighborhood probability
+
+            // Loop over the neighborhoods
             for(j=0; j<conf_info.get_n_nbrhd(); j++) {
 
-               smooth_dp = smooth_field(prob_dp, InterpMthd_UW_Mean,
+               nbrhd_dp = smooth_field(prob_dp, InterpMthd_UW_Mean,
                               conf_info.nbrhd_prob.width[j],
-                              conf_info.nbrhd_prob.shape,
-                              grid.wrap_lon(),
+                              conf_info.nbrhd_prob.shape, grid.wrap_lon(),
                               conf_info.nbrhd_prob.vld_thresh, info);
 
-               for(k=0; k<count_na.n(); k++) {
-                  ens_prob[k] = (float) smooth_dp.buf()[k];
-               }
-
-               // Output variable name
+               // Write neighborhood ensemble probability
                snprintf(type_str, sizeof(type_str), "ENS_NEP_%s_%s%i",
-                        conf_info.ens_ta[i_ens][i].get_abbr_str().contents().c_str(),
+                        conf_info.ens_ta[i_var][i].get_abbr_str().contents().c_str(),
                         interpmthd_to_string(InterpMthd_Nbrhd).c_str(),
                         conf_info.nbrhd_prob.width[j]*conf_info.nbrhd_prob.width[j]);
-
-               write_ens_var_float(i_ens, ens_prob, dp,
-                                   type_str,
-                                   "Neighborhood Ensemble Probability");
-
+               write_ens_data_plane(i_var, nbrhd_dp, ens_dp, type_str,
+                                    "Neighborhood Ensemble Probability");
             } // end for j
          } // end if do_nep
       } // end for i
@@ -2523,10 +2506,10 @@ void write_ens_nc(int i_ens, DataPlane &dp) {
    // Add the neighborhood maximum ensemble probabilities, if requested
    if(conf_info.nc_info.do_nmep) {
 
-      prob_dp.set_size(dp.nx(), dp.ny());
+      prob_dp.set_size(grid.nx(), grid.ny());
 
       // Loop through each threshold
-      for(i=0; i<conf_info.ens_ta[i_ens].n(); i++) {
+      for(i=0; i<conf_info.ens_ta[i_var].n(); i++) {
 
          // Loop through each neigbhorhood size
          for(j=0; j<conf_info.get_n_nbrhd(); j++) {
@@ -2535,44 +2518,31 @@ void write_ens_nc(int i_ens, DataPlane &dp) {
             prob_dp.erase();
 
             // Compute the neighborhood maximum ensemble probability
-            for(k=0; k<count_na.n(); k++) {
-
-               // Check for too much missing data
-               if((double) (count_na[k]/n_ens_vld[i_ens]) < t) {
-                  prob_dp.buf()[k] = bad_data_double;
-               }
-               else {
-                  prob_dp.buf()[k] = (double) (thresh_nbrhd_count_na[i][j][k]/count_na[k]);
-               }
-            } // end for k
+            for(k=0; k<cnt_na.n(); k++) {
+               prob_dp.buf()[k] = ((double) (cnt_na[k]/n_ens_vld[i_var]) < t ?
+                                   bad_data_double :
+                                   (double) (thresh_nbrhd_cnt_na[i][j][k]/cnt_na[k]));
+            }
 
             // Apply requested NMEP smoothers
             for(k=0; k<conf_info.nmep_smooth.n_interp; k++) {
 
-               smooth_dp = smooth_field(prob_dp,
-                              string_to_interpmthd(conf_info.nmep_smooth.method[k].c_str()),
-                              conf_info.nmep_smooth.width[k],
-                              conf_info.nmep_smooth.shape,
-                              grid.wrap_lon(),
-                              conf_info.nmep_smooth.vld_thresh,
-                              conf_info.nmep_smooth.gaussian);
-
-               for(l=0; l<count_na.n(); l++) {
-                  ens_prob[l] = (float) smooth_dp.buf()[l];
-               }
+               nbrhd_dp = smooth_field(prob_dp,
+                             string_to_interpmthd(conf_info.nmep_smooth.method[k].c_str()),
+                             conf_info.nmep_smooth.width[k],
+                             conf_info.nmep_smooth.shape, grid.wrap_lon(),
+                             conf_info.nmep_smooth.vld_thresh,
+                             conf_info.nmep_smooth.gaussian);
 
                // Output variable name
                snprintf(type_str, sizeof(type_str), "ENS_NMEP_%s_%s%i_%s%i",
-                        conf_info.ens_ta[i_ens][i].get_abbr_str().contents().c_str(),
+                        conf_info.ens_ta[i_var][i].get_abbr_str().contents().c_str(),
                         interpmthd_to_string(InterpMthd_Nbrhd).c_str(),
                         conf_info.nbrhd_prob.width[j]*conf_info.nbrhd_prob.width[j],
                         conf_info.nmep_smooth.method[k].c_str(),
                         conf_info.nmep_smooth.width[k]*conf_info.nmep_smooth.width[k]);
-
-               write_ens_var_float(i_ens, ens_prob, dp,
-                                   type_str,
-                                   "Neighborhood Maximum Ensemble Probability");
-
+               write_ens_data_plane(i_var, nbrhd_dp, ens_dp, type_str,
+                                    "Neighborhood Maximum Ensemble Probability");
             } // end for k
          } // end for j
       } // end for i
@@ -2587,27 +2557,26 @@ void write_ens_nc(int i_ens, DataPlane &dp) {
    if(ens_max)   { delete [] ens_max;   ens_max   = (float *) 0; }
    if(ens_range) { delete [] ens_range; ens_range = (float *) 0; }
    if(ens_vld)   { delete [] ens_vld;   ens_vld   = (int   *) 0; }
-   if(ens_prob)  { delete [] ens_prob;  ens_prob  = (float *) 0; }
 
    return;
 }
 
 ////////////////////////////////////////////////////////////////////////
 
-void write_ens_var_float(int i_ens, float *ens_data, DataPlane &dp,
+void write_ens_var_float(int i_var, float *ens_data, const DataPlane &dp,
                          const char *type_str,
                          const char *long_name_str) {
    NcVar ens_var;
    ConcatString ens_var_name, var_str, name_str, cs;
 
    // Append nc_pairs_var_str config file entry
-   cs = conf_info.ens_var_str[i_ens];
+   cs = conf_info.ens_var_str[i_var];
    if(cs.length() > 0) var_str << "_" << cs;
 
    // Construct the variable name
    ens_var_name << cs_erase
-                << conf_info.ens_info[i_ens]->name_attr() << "_"
-                << conf_info.ens_info[i_ens]->level_attr()
+                << conf_info.ens_info[i_var]->name_attr() << "_"
+                << conf_info.ens_info[i_var]->level_attr()
                 << var_str << "_" << type_str;
 
    // Skip variable names that have already been written
@@ -2625,16 +2594,16 @@ void write_ens_var_float(int i_ens, float *ens_data, DataPlane &dp,
    //
    if(strcmp(type_str, "ENS_MEAN") == 0) {
       name_str << cs_erase
-               << conf_info.ens_info[i_ens]->name_attr();
+               << conf_info.ens_info[i_var]->name_attr();
    }
    else {
       name_str << cs_erase
-               << conf_info.ens_info[i_ens]->name_attr() << "_"
+               << conf_info.ens_info[i_var]->name_attr() << "_"
                << type_str;
    }
 
    // Add the variable attributes
-   add_var_att_local(conf_info.ens_info[i_ens], &ens_var, false, dp,
+   add_var_att_local(conf_info.ens_info[i_var], &ens_var, false, dp,
                      name_str.c_str(), long_name_str);
 
    // Write the data
@@ -2650,20 +2619,20 @@ void write_ens_var_float(int i_ens, float *ens_data, DataPlane &dp,
 
 ////////////////////////////////////////////////////////////////////////
 
-void write_ens_var_int(int i_ens, int *ens_data, DataPlane &dp,
+void write_ens_var_int(int i_var, int *ens_data, const DataPlane &dp,
                        const char *type_str,
                        const char *long_name_str) {
    NcVar ens_var;
    ConcatString ens_var_name, var_str, name_str, cs;
 
    // Append nc_pairs_var_str config file entry
-   cs = conf_info.ens_var_str[i_ens];
+   cs = conf_info.ens_var_str[i_var];
    if(cs.length() > 0) var_str << "_" << cs;
 
    // Construct the variable name
    ens_var_name << cs_erase
-                << conf_info.ens_info[i_ens]->name_attr() << "_"
-                << conf_info.ens_info[i_ens]->level_attr()
+                << conf_info.ens_info[i_var]->name_attr() << "_"
+                << conf_info.ens_info[i_var]->level_attr()
                 << var_str << "_" << type_str;
 
    // Skip variable names that have already been written
@@ -2678,11 +2647,11 @@ void write_ens_var_int(int i_ens, int *ens_data, DataPlane &dp,
 
    // Construct the variable name attribute
    name_str << cs_erase
-            << conf_info.ens_info[i_ens]->name_attr() << "_"
+            << conf_info.ens_info[i_var]->name_attr() << "_"
             << type_str;
 
    // Add the variable attributes
-   add_var_att_local(conf_info.ens_info[i_ens], &ens_var, true, dp,
+   add_var_att_local(conf_info.ens_info[i_var], &ens_var, true, dp,
                      name_str.c_str(), long_name_str);
 
    // Write the data
@@ -2698,6 +2667,26 @@ void write_ens_var_int(int i_ens, int *ens_data, DataPlane &dp,
 
 ////////////////////////////////////////////////////////////////////////
 
+void write_ens_data_plane(int i_var, const DataPlane &ens_dp, const DataPlane &dp,
+                          const char *type_str,  const char *long_name_str) {
+
+   // Allocate memory for this data
+   float *ens_data = new float [nxy];
+
+   // Store the data in an array of floats
+   for(int i=0; i<nxy; i++) ens_data[i] = ens_dp.data()[i];
+
+   // Write the output
+   write_ens_var_float(i_var, ens_data, dp, type_str, long_name_str);
+
+   // Cleanup
+   if(ens_data) { delete [] ens_data; ens_data = (float *) 0; }
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void write_orank_nc(PairDataEnsemble &pd, DataPlane &dp,
                     int i_vx, int i_interp, int i_mask) {
    int i, n;
@@ -2709,13 +2698,13 @@ void write_orank_nc(PairDataEnsemble &pd, DataPlane &dp,
    int   *ens_vld  = (int *)   0;
 
    // Allocate memory for storing ensemble data
-   obs_v    = new float [grid.nx()*grid.ny()];
-   obs_rank = new int   [grid.nx()*grid.ny()];
-   obs_pit  = new float [grid.nx()*grid.ny()];
-   ens_vld  = new int   [grid.nx()*grid.ny()];
+   obs_v    = new float [nxy];
+   obs_rank = new int   [nxy];
+   obs_pit  = new float [nxy];
+   ens_vld  = new int   [nxy];
 
    // Initialize
-   for(i=0; i<grid.nx()*grid.ny(); i++) {
+   for(i=0; i<nxy; i++) {
       obs_v[i]    = bad_data_float;
       obs_rank[i] = bad_data_int;
       obs_pit[i]  = bad_data_float;
@@ -2895,8 +2884,10 @@ void write_orank_var_int(int i_vx, int i_interp, int i_mask,
 
 ////////////////////////////////////////////////////////////////////////
 
-void add_var_att_local(VarInfo *info, NcVar *nc_var, bool is_int, DataPlane &dp,
-                       const char *name_str, const char *long_name_str) {
+void add_var_att_local(VarInfo *info, NcVar *nc_var, bool is_int,
+                       const DataPlane &dp,
+                       const char *name_str,
+                       const char *long_name_str) {
    ConcatString att_str;
 
    // Construct the long name
@@ -2966,23 +2957,23 @@ void clean_up() {
    }
 
    // Deallocate threshold count arrays
-   if(thresh_count_na) {
+   if(thresh_cnt_na) {
       for(i=0; i<conf_info.get_max_n_thresh(); i++) {
-         thresh_count_na[i].clear();
+         thresh_cnt_na[i].clear();
       }
-      delete [] thresh_count_na;
-      thresh_count_na = (NumArray *) 0;
+      delete [] thresh_cnt_na;
+      thresh_cnt_na = (NumArray *) 0;
    }
-   if(thresh_nbrhd_count_na) {
+   if(thresh_nbrhd_cnt_na) {
       for(i=0; i<conf_info.get_max_n_thresh(); i++) {
          for(j=0; j<conf_info.get_n_nbrhd(); j++) {
-            thresh_nbrhd_count_na[i][j].clear();
+            thresh_nbrhd_cnt_na[i][j].clear();
          }
-         delete [] thresh_nbrhd_count_na[i];
-         thresh_nbrhd_count_na[i] = (NumArray *) 0;
+         delete [] thresh_nbrhd_cnt_na[i];
+         thresh_nbrhd_cnt_na[i] = (NumArray *) 0;
       }
-      delete [] thresh_nbrhd_count_na;
-      thresh_nbrhd_count_na = (NumArray **) 0;
+      delete [] thresh_nbrhd_cnt_na;
+      thresh_nbrhd_cnt_na = (NumArray **) 0;
    }
 
    return;
@@ -3001,6 +2992,7 @@ void usage() {
         << "\t[-grid_obs file]\n"
         << "\t[-point_obs file]\n"
         << "\t[-ens_mean file]\n"
+        << "\t[-ctrl file]\n"
         << "\t[-obs_valid_beg time]\n"
         << "\t[-obs_valid_end time]\n"
         << "\t[-outdir path]\n"
@@ -3027,6 +3019,9 @@ void usage() {
         << "\t\t\"-ens_mean file\" specifies an ensemble mean model data file "
         << "(optional).\n"
 
+        << "\t\t\"-ctrl file\" specifies the gridded ensemble control data file "
+        << "included in the mean but excluded from the spread (optional).\n"
+
         << "\t\t\"-obs_valid_beg time\" in YYYYMMDD[_HH[MMSS]] sets the "
         << "beginning of the matching time window (optional).\n"
 
@@ -3073,6 +3068,12 @@ void set_ens_mean(const StringArray & a)
 
 ////////////////////////////////////////////////////////////////////////
 
+void set_ctrl_file(const StringArray & a) {
+   ctrl_file = a[0];
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void set_obs_valid_beg(const StringArray & a)
 {
    obs_valid_beg_ut = timestring_to_unix(a[0].c_str());
diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.h b/met/src/tools/core/ensemble_stat/ensemble_stat.h
index c3b9feec37..784eecad0e 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.h
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.h
@@ -105,6 +105,8 @@ static bool         ens_mean_flag; // Flag for ensemble mean processing
 static ConcatString ens_mean_user; // User-specified ensemble mean data file
 static ConcatString ens_mean_file; // Computed ensemble mean output file
 
+static ConcatString ctrl_file;     // Control member
+
 // Input Observation files
 static StringArray  grid_obs_file_list;
 static int          grid_obs_flag = 0;
@@ -121,11 +123,11 @@ static ConcatString         config_file;
 static ConcatString         out_file;
 
 // Optional arguments
-static unixtime     ens_valid_ut        = (unixtime) 0;
-static NumArray     ens_lead_na;
+static unixtime ens_valid_ut = (unixtime) 0;
+static NumArray ens_lead_na;
 
-static unixtime     obs_valid_beg_ut    = (unixtime) 0;
-static unixtime     obs_valid_end_ut    = (unixtime) 0;
+static unixtime obs_valid_beg_ut = (unixtime) 0;
+static unixtime obs_valid_end_ut = (unixtime) 0;
 
 static ConcatString out_dir;
 
@@ -179,9 +181,10 @@ static Met2dDataFileFactory mtddf_factory;
 static StatHdrColumns shc;
 
 // Arrays to store running sums and counts
-static NumArray count_na, min_na, max_na, sum_na, sum_sq_na;
-static NumArray *thresh_count_na = (NumArray *) 0; // [n_thresh]
-static NumArray **thresh_nbrhd_count_na = (NumArray **) 0; // [n_thresh][n_nbrhd]
+static NumArray cnt_na, min_na, max_na, sum_na;
+static NumArray stdev_cnt_na, stdev_sum_na, stdev_ssq_na;
+static NumArray *thresh_cnt_na = (NumArray *) 0; // [n_thresh]
+static NumArray **thresh_nbrhd_cnt_na = (NumArray **) 0; // [n_thresh][n_nbrhd]
 
 ////////////////////////////////////////////////////////////////////////
 

From b39d9de96ca1cbb313d6962dfd189494961736e8 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 8 Oct 2021 11:49:45 -0600
Subject: [PATCH 03/19] Per #1905, update the users guide with the -ctrl
 command line option.

---
 met/docs/Users_Guide/ensemble-stat.rst | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst
index a39424d289..cf2cbd0012 100644
--- a/met/docs/Users_Guide/ensemble-stat.rst
+++ b/met/docs/Users_Guide/ensemble-stat.rst
@@ -69,6 +69,7 @@ The usage statement for the Ensemble Stat tool is shown below:
          [-grid_obs file]
          [-point_obs file]
          [-ens_mean file]
+         [-ctrl file]
          [-obs_valid_beg time]
          [-obs_valid_end time]
          [-outdir path]
@@ -92,23 +93,23 @@ Optional arguments for ensemble_stat
 
 4. To produce ensemble statistics using gridded observations, use the **-grid_obs file** option to specify a gridded observation file. This option may be used multiple times if your observations are in several files.
 
+5. To produce ensemble statistics using point observations, use the **-point_obs file** option to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files.
 
-5. To produce ensemble statistics using point observations, use the **-point_obs file** to specify a NetCDF point observation file. This option may be used multiple times if your observations are in several files.
+6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** option specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** option from earlier versions of MET.
 
+7. The **-ctrl file** option specifies and ensemble control member data file. The control member is included in the computation of the ensemble mean but excluded from the spread.
 
-6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** from earlier versions of MET.
+8. To filter point observations by time, use **-obs_valid_beg time** in YYYYMMDD[_HH[MMSS]] format to set the beginning of the matching observation time window.
 
-7. To filter point observations by time, use **-obs_valid_beg time** in YYYYMMDD[_HH[MMSS]] format to set the beginning of the matching observation time window.
+9. As above, use **-obs_valid_end time** in YYYYMMDD[_HH[MMSS]] format to set the end of the matching observation time window.
 
-8. As above, use **-obs_valid_end time** in YYYYMMDD[_HH[MMSS]] format to set the end of the matching observation time window.
+10. Specify the **-outdir path** option to override the default output directory (./).
 
-9. Specify the **-outdir path** option to override the default output directory (./).
+11. The **-log** file outputs log messages to the specified file.
 
-10. The **-log** file outputs log messages to the specified file.
+12. The **-v level** option indicates the desired level of verbosity. The value of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity will increase the amount of logging.
 
-11. The **-v level** option indicates the desired level of verbosity. The value of "level" will override the default setting of 2. Setting the verbosity to 0 will make the tool run with no log messages, while increasing the verbosity will increase the amount of logging.
-
-12. The **-compress level** option indicates the desired level of compression (deflate level) for NetCDF variables. The valid level is between 0 and 9. The value of "level" will override the default setting of 0 from the configuration file or the environment variable MET_NC_COMPRESS. Setting the compression level to 0 will make no compression for the NetCDF output. Lower number is for fast compression and higher number is for better compression.
+13. The **-compress level** option indicates the desired level of compression (deflate level) for NetCDF variables. The valid level is between 0 and 9. The value of "level" will override the default setting of 0 from the configuration file or the environment variable MET_NC_COMPRESS. Setting the compression level to 0 will make no compression for the NetCDF output. Lower number is for fast compression and higher number is for better compression.
 
 An example of the ensemble_stat calling sequence is shown below:
 

From ae054b1aa4c97dadcd3f20758e5c75b6ce754266 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Mon, 11 Oct 2021 10:18:48 -0600
Subject: [PATCH 04/19] Not actually related to #1905, but clarifying the
 pcp_combine usage statement per dtcenter/METplus#1201 discussion.

---
 met/docs/Users_Guide/reformat_grid.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/met/docs/Users_Guide/reformat_grid.rst b/met/docs/Users_Guide/reformat_grid.rst
index 191ef4e1a4..0c186d524e 100644
--- a/met/docs/Users_Guide/reformat_grid.rst
+++ b/met/docs/Users_Guide/reformat_grid.rst
@@ -74,7 +74,7 @@ Required arguments for the pcp_combine
 Optional arguments for pcp_combine
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. This option may be used multiple times and output will be created for each.
+3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the name and level of the requested data and be enclosed in single quotes (e.g. -field 'name="precip06h"; level="(*,*)";' for reading a NetCDF variable). It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options.
 
 4. The **-name list** option is a comma-separated list of output variable names which override the default choices. If specified, the number of names must match the number of variables to be written to the output file.
 

From f2339a3669977cfb85658df7a5b138edabfb6b97 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Mon, 11 Oct 2021 10:29:40 -0600
Subject: [PATCH 05/19] Avoid inline stars in pcp_combine example that don't
 render well.

---
 met/docs/Users_Guide/reformat_grid.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/met/docs/Users_Guide/reformat_grid.rst b/met/docs/Users_Guide/reformat_grid.rst
index 0c186d524e..d4f7190e70 100644
--- a/met/docs/Users_Guide/reformat_grid.rst
+++ b/met/docs/Users_Guide/reformat_grid.rst
@@ -74,7 +74,7 @@ Required arguments for the pcp_combine
 Optional arguments for pcp_combine
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the name and level of the requested data and be enclosed in single quotes (e.g. -field 'name="precip06h"; level="(*,*)";' for reading a NetCDF variable). It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options.
+3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the name and level of the requested data and be enclosed in single quotes (e.g. -field 'name="ACPCP"; level="A6";' for 6-hourly accumulated convective precipitation from a GRIB file). It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options.
 
 4. The **-name list** option is a comma-separated list of output variable names which override the default choices. If specified, the number of names must match the number of variables to be written to the output file.
 

From 52935a9661f3a06fa48c8d552039a2229025c558 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Mon, 11 Oct 2021 10:40:16 -0600
Subject: [PATCH 06/19] Per #1905, tweak the docs.

---
 met/docs/Users_Guide/reformat_grid.rst | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/met/docs/Users_Guide/reformat_grid.rst b/met/docs/Users_Guide/reformat_grid.rst
index d4f7190e70..143cdcbeb4 100644
--- a/met/docs/Users_Guide/reformat_grid.rst
+++ b/met/docs/Users_Guide/reformat_grid.rst
@@ -74,7 +74,9 @@ Required arguments for the pcp_combine
 Optional arguments for pcp_combine
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the name and level of the requested data and be enclosed in single quotes (e.g. -field 'name="ACPCP"; level="A6";' for 6-hourly accumulated convective precipitation from a GRIB file). It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options.
+3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each.
+
+In general, the field string should include the name and level of the requested data and be enclosed in single quotes. It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options. For example, use **-field ‘name=”ACPCP”; level=”A6”; convert(x)=x/25.4;’ to read 6-hourly accumulated convective precipitation from a GRIB file and convert from millimeters to inches. 
 
 4. The **-name list** option is a comma-separated list of output variable names which override the default choices. If specified, the number of names must match the number of variables to be written to the output file.
 

From e89f122334fc992018c757daa9f0bf48f2f99705 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Mon, 11 Oct 2021 10:46:51 -0600
Subject: [PATCH 07/19] Per #1905, more tweaks.

---
 met/docs/Users_Guide/reformat_grid.rst | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/met/docs/Users_Guide/reformat_grid.rst b/met/docs/Users_Guide/reformat_grid.rst
index 143cdcbeb4..3f67a377a7 100644
--- a/met/docs/Users_Guide/reformat_grid.rst
+++ b/met/docs/Users_Guide/reformat_grid.rst
@@ -74,9 +74,7 @@ Required arguments for the pcp_combine
 Optional arguments for pcp_combine
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each.
-
-In general, the field string should include the name and level of the requested data and be enclosed in single quotes. It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options. For example, use **-field ‘name=”ACPCP”; level=”A6”; convert(x)=x/25.4;’ to read 6-hourly accumulated convective precipitation from a GRIB file and convert from millimeters to inches. 
+3. The **-field string** option defines the data to be extracted from the input files. Use this option when processing fields other than **APCP** or non-GRIB files. It can be used multiple times and output will be created for each. In general, the field string should include the **name** and **level** of the requested data and be enclosed in single quotes. It is processed as an inline configuration file and may also include data filtering, censoring, and conversion options. For example, use **-field ‘name=”ACPCP”; level=”A6”; convert(x)=x/25.4;’** to read 6-hourly accumulated convective precipitation from a GRIB file and convert from millimeters to inches. 
 
 4. The **-name list** option is a comma-separated list of output variable names which override the default choices. If specified, the number of names must match the number of variables to be written to the output file.
 

From 2a63c70353035a2808479a6146e9b1bfa0b84d1e Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Wed, 27 Oct 2021 10:15:47 -0600
Subject: [PATCH 08/19] Per #1905, just aligning code.

---
 .../tools/core/ensemble_stat/ensemble_stat.cc | 21 ++++++++-----------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
index 77420adca5..fc39e16776 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -219,21 +219,18 @@ void process_command_line(int argc, char **argv) {
    cline.set_usage(usage);
 
    //
-   // add the options function calls
+   // add the option function calls
+   // quietly support deprecated -ssvar_mean option
    //
-   cline.add(set_grid_obs,      "-grid_obs", 1);
-   cline.add(set_point_obs,     "-point_obs", 1);
-   cline.add(set_ens_mean,      "-ens_mean", 1);
-   cline.add(set_ctrl_file,     "-ctrl", 1);
+   cline.add(set_grid_obs,      "-grid_obs",      1);
+   cline.add(set_point_obs,     "-point_obs",     1);
+   cline.add(set_ens_mean,      "-ens_mean",      1);
+   cline.add(set_ctrl_file,     "-ctrl",          1);
    cline.add(set_obs_valid_beg, "-obs_valid_beg", 1);
    cline.add(set_obs_valid_end, "-obs_valid_end", 1);
-   cline.add(set_outdir,        "-outdir", 1);
-   cline.add(set_compress,      "-compress",  1);
-
-   //
-   // quietly support deprecated -ssvar_mean option
-   //
-   cline.add(set_ens_mean,      "-ssvar_mean", 1);
+   cline.add(set_outdir,        "-outdir",        1);
+   cline.add(set_compress,      "-compress",      1);
+   cline.add(set_ens_mean,      "-ssvar_mean",    1);
 
    //
    // parse the command line

From 519daad5680168cef23836c3e735728cfa50968c Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Wed, 27 Oct 2021 10:27:59 -0600
Subject: [PATCH 09/19] Per #1905, add a call to ensemble_stat using the -ctrl
 option. I realize that if a control member is specified, the number of
 ensemble members should be increased by 1.

---
 test/xml/unit_ensemble_stat.xml | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/test/xml/unit_ensemble_stat.xml b/test/xml/unit_ensemble_stat.xml
index ab7cfdef42..105d19c3d4 100644
--- a/test/xml/unit_ensemble_stat.xml
+++ b/test/xml/unit_ensemble_stat.xml
@@ -122,6 +122,34 @@
     </output>
   </test>
 
+  <test name="ensemble_stat_MASK_SID_CTRL">
+    <exec>&MET_BIN;/ensemble_stat</exec>
+    <env>
+      <pair><name>CENSOR_THRESH</name>  <value></value></pair>
+      <pair><name>CENSOR_VAL</name>     <value></value></pair>
+      <pair><name>SKIP_CONST</name>     <value>FALSE</value></pair>
+      <pair><name>OUTPUT_PREFIX</name>  <value>MASK_SID_CTRL</value></pair>
+      <pair><name>CONFIG_DIR</name>     <value>&CONFIG_DIR;</value></pair>
+    </env>
+    <param> \
+      5 \
+      &DATA_DIR_MODEL;/grib1/arw-fer-gep5/arw-fer-gep5_2012040912_F024.grib \
+      &DATA_DIR_MODEL;/grib1/arw-sch-gep2/arw-sch-gep2_2012040912_F024.grib \
+      &DATA_DIR_MODEL;/grib1/arw-sch-gep6/arw-sch-gep6_2012040912_F024.grib \
+      &DATA_DIR_MODEL;/grib1/arw-tom-gep3/arw-tom-gep3_2012040912_F024.grib \
+      &DATA_DIR_MODEL;/grib1/arw-tom-gep7/arw-tom-gep7_2012040912_F024.grib \
+      &CONFIG_DIR;/EnsembleStatConfig_MASK_SID \
+      -ctrl &DATA_DIR_MODEL;/grib1/arw-fer-gep1/arw-fer-gep1_2012040912_F024.grib \
+      -point_obs &OUTPUT_DIR;/ascii2nc/gauge_2012041012_24hr.nc \
+      -outdir &OUTPUT_DIR;/ensemble_stat -v 1
+    </param>
+    <output>
+      <stat>&OUTPUT_DIR;/ensemble_stat/ensemble_stat_MASK_SID_CTRL_20120410_120000V.stat</stat>
+      <stat>&OUTPUT_DIR;/ensemble_stat/ensemble_stat_MASK_SID_CTRL_20120410_120000V_orank.txt</stat>
+      <grid_nc>&OUTPUT_DIR;/ensemble_stat/ensemble_stat_MASK_SID_CTRL_20120410_120000V_ens.nc</grid_nc>
+    </output>
+  </test>
+
   <test name="ensemble_stat_MASK_SID_CENSOR">
     <exec>&MET_BIN;/ensemble_stat</exec>
     <env>

From 1e25123cb9bf9dad5c1bb24307031d7c5b623e1a Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Wed, 27 Oct 2021 15:38:57 -0600
Subject: [PATCH 10/19] Per #1905, enable STDEV output in the NetCDF file since
 it's computation differs when -ctrl is used.

---
 test/config/EnsembleStatConfig_MASK_SID | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/config/EnsembleStatConfig_MASK_SID b/test/config/EnsembleStatConfig_MASK_SID
index d33e7d9bd9..2ad34e2f27 100644
--- a/test/config/EnsembleStatConfig_MASK_SID
+++ b/test/config/EnsembleStatConfig_MASK_SID
@@ -206,7 +206,7 @@ output_flag = {
 //
 ensemble_flag = {
    mean      = TRUE;
-   stdev     = FALSE;
+   stdev     = TRUE;
    minus     = FALSE;
    plus      = FALSE;
    min       = TRUE;

From e92e04b4ae7a7f4510b80ef3433fe72fec89dd76 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 16:31:44 -0600
Subject: [PATCH 11/19] Per #1905, add NumArray::mean(), variance(), and
 stdev() with options to skip. This supports excluding the ensemble control
 member from the spread.

---
 met/src/basic/vx_util/num_array.cc | 53 ++++++++++++++++++++++++++++++
 met/src/basic/vx_util/num_array.h  |  4 +++
 2 files changed, 57 insertions(+)

diff --git a/met/src/basic/vx_util/num_array.cc b/met/src/basic/vx_util/num_array.cc
index b681b901ee..446bf29034 100644
--- a/met/src/basic/vx_util/num_array.cc
+++ b/met/src/basic/vx_util/num_array.cc
@@ -1157,6 +1157,59 @@ double NumArray::wmean_fisher(const NumArray &wgt) const
 
 }
 
+
+////////////////////////////////////////////////////////////////////////
+
+
+double NumArray::variance(int skip_index) const
+
+{
+
+   if(n() == 0)  return ( bad_data_double );
+
+   int j, count;
+   double s, s_sq, var;
+
+   s = s_sq = 0.0;
+   count = 0;
+
+   for(j=0; j<n(); j++) {
+      if(is_bad_data(e[j]) || j == skip_index) continue;
+      s    += e[j];
+      s_sq += e[j]*e[j];
+      count++;
+   }
+
+   // Check for slightly negative precision error
+   if(count > 1) {
+      var = (s_sq - s*s/(double) count)/((double) (count - 1));
+      if(is_eq(var, 0.0)) var = 0.0;
+   }
+   else {
+      var = bad_data_double;
+   }
+
+   return(var);
+
+}
+
+
+////////////////////////////////////////////////////////////////////////
+
+
+double NumArray::stdev(int skip_index) const
+
+{
+
+   double v = variance(skip_index);
+
+   if ( !is_bad_data(v) )  v = sqrt(v);
+
+   return(v);
+
+}
+
+
 ////////////////////////////////////////////////////////////////////////
 
 //
diff --git a/met/src/basic/vx_util/num_array.h b/met/src/basic/vx_util/num_array.h
index 207b512acd..623f7e39c8 100644
--- a/met/src/basic/vx_util/num_array.h
+++ b/met/src/basic/vx_util/num_array.h
@@ -21,6 +21,7 @@
 #include <iostream>
 
 #include "concat_string.h"
+#include "is_bad_data.h"
 
 
 ////////////////////////////////////////////////////////////////////////
@@ -107,6 +108,9 @@ class NumArray {
       double mean_sqrt() const;
       double mean_fisher() const;
 
+      double variance(int skip_index = bad_data_int) const;
+      double stdev(int skip_index = bad_data_int) const;
+
       double wmean(const NumArray &) const;
       double wmean_sqrt(const NumArray &) const;
       double wmean_fisher(const NumArray &) const;

From a8e6ca8270f6385dd98904a291c6d75f1977956d Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 16:36:03 -0600
Subject: [PATCH 12/19] Per #1905, keep track of which ensemble member is the
 control and exclude it from the spread and variance computations.

---
 .../vx_statistics/pair_data_ensemble.cc       | 57 ++++++++++++++-----
 .../vx_statistics/pair_data_ensemble.h        |  4 ++
 2 files changed, 46 insertions(+), 15 deletions(-)

diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.cc b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
index 5782cf6a0d..77b66536b3 100644
--- a/met/src/libcode/vx_statistics/pair_data_ensemble.cc
+++ b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
@@ -108,6 +108,7 @@ void PairDataEnsemble::clear() {
 
    n_ens = 0;
    n_pair = 0;
+   ctrl_index = bad_data_int;
    skip_const = false;
    skip_ba.clear();
 
@@ -121,6 +122,7 @@ void PairDataEnsemble::clear() {
 
    esum_na.clear();
    esumsq_na.clear();
+   esumn_na.clear();
 
    mn_na.clear();
    mn_oerr_na.clear();
@@ -170,6 +172,7 @@ void PairDataEnsemble::extend(int n) {
    var_plus_oerr_na.extend   (n);
    esum_na.extend            (n);
    esumsq_na.extend          (n);
+   esumn_na.extend           (n);
    mn_na.extend              (n);
    mn_oerr_na.extend         (n);
 
@@ -222,6 +225,7 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) {
    ign_na         = pd.ign_na;
    pit_na         = pd.pit_na;
    n_pair         = pd.n_pair;
+   ctrl_index     = pd.ctrl_index;
    skip_const     = pd.skip_const;
    skip_ba        = pd.skip_ba;
 
@@ -231,6 +235,7 @@ void PairDataEnsemble::assign(const PairDataEnsemble &pd) {
 
    esum_na        = pd.esum_na;
    esumsq_na      = pd.esumsq_na;
+   esumn_na       = pd.esumn_na;
 
    mn_na          = pd.mn_na;
    mn_oerr_na     = pd.mn_oerr_na;
@@ -290,12 +295,14 @@ void PairDataEnsemble::add_ens_var_sums(int i_obs, double v) {
    if(i_obs >= esum_na.n()) {
       esum_na.add(0.0);
       esumsq_na.add(0.0);
+      esumn_na.add(0.0);
    }
 
    // Track sums of the raw ensemble member values
    if(!is_bad_data(v)) {
       esum_na.inc(i_obs, v);
       esumsq_na.inc(i_obs, v*v);
+      esumn_na.inc(i_obs, 1);
    }
 
    return;
@@ -416,17 +423,17 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
       else {
 
          // Compute the variance of the unperturbed ensemble members
-         var_unperturbed = compute_variance(esum_na[i], esumsq_na[i], v_na[i]);
+         var_unperturbed = compute_variance(esum_na[i], esumsq_na[i], esumn_na[i]);
          var_na.add(var_unperturbed);
 
          // Process the observation error information.
          ObsErrorEntry * e = (has_obs_error() ? obs_error_entry[i] : 0);
          if(e) {
 
-            // Compute perturbed ensemble mean and variance 
-            cur_ens.compute_mean_variance(mean, var_perturbed);
-            mn_oerr_na.add(mean);
-            var_oerr_na.add(var_perturbed);
+            // Compute perturbed ensemble mean and variance
+            // Exclude the control member from the variance
+            mn_oerr_na.add(cur_ens.mean());
+            var_oerr_na.add(cur_ens.variance(ctrl_index)); // JHG
 
             // Compute the variance plus observation error variance.
             var_plus_oerr_na.add(var_unperturbed +
@@ -467,8 +474,12 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
          crps_emp_na.add(compute_crps_emp(o_na[i], cur_ens));
          crpscl_emp_na.add(compute_crps_emp(o_na[i], cur_clm));
 
+         // Ensemble mean and standard deviation
+         // Exclude the control member from the standard deviation
+         mean = cur_ens.mean();
+         stdev = cur_ens.stdev(ctrl_index); // JHG
+
          // Store Gaussian CRPS stats
-         cur_ens.compute_mean_stdev(mean, stdev);
          crps_gaus_na.add(compute_crps_gaus(o_na[i], mean, stdev));
          crpscl_gaus_na.add(compute_crps_gaus(o_na[i], cmn_na[i], csd_na[i]));
          ign_na.add(compute_ens_ign(o_na[i], mean, stdev));
@@ -613,7 +624,7 @@ struct ssvar_bin_comp {
 
 void PairDataEnsemble::compute_ssvar() {
    int i, j;
-   double mn, var;
+   double var;
    ssvar_bin_map bins;
    NumArray cur;
 
@@ -641,24 +652,22 @@ void PairDataEnsemble::compute_ssvar() {
    for(i=0; i<o_na.n(); i++) {
 
       // Check if point should be skipped
-      if(skip_ba[i]) continue;
+      // Exclude the control member from the variance
+      if(skip_ba[i] || i == ctrl_index) continue;
 
       // Store ensemble values for the current observation
       for(j=0, cur.erase(); j<n_ens; j++) cur.add(e_na[j][i]);
 
-      // Compute standard deviation
-      cur.compute_mean_variance(mn, var);
-
       // Build a variance point
       ens_ssvar_pt pt;
-      pt.var = var;
+      pt.var = cur.variance();
       pt.f   = mn_na[i];
       pt.o   = o_na[i];
       pt.w   = wgt_na[i];
 
       // Determine the bin for the current point and add it to the list
       // Bins are defined starting at 0 and are left-closed, right-open
-      j = floor(var/ssvar_bin_size);
+      j = floor(pt.var/ssvar_bin_size);
       string ssvar_min = str_format("%.5e", j*ssvar_bin_size).contents();
       if( !bins.count(ssvar_min) ){
          ssvar_pt_list pts;
@@ -801,7 +810,7 @@ PairDataEnsemble PairDataEnsemble::subset_pairs_obs_thresh(const SingleThresh &o
       //
       // Exclude from subset:
       //   sid_sa, lat_na, lon_na, x_na, y_na, vld_ta, lvl_ta, elv_ta,
-      //   o_qc_sa, esum_na, esumsq_na
+      //   o_qc_sa, esum_na, esumsq_na, esumn_na
 
       pd.wgt_na.add(wgt_na[i]);
       pd.o_na.add(o_na[i]);
@@ -1667,7 +1676,10 @@ void VxPairDataEnsemble::add_ens(int member, bool mn, Grid &gr) {
                else {
 
                   // Track unperturbed ensemble variance sums
-                  pd[i][j][k].add_ens_var_sums(l, fcst_v);
+                  // Exclude the control member from the variance
+                  if(member != pd[i][j][k].ctrl_index) {
+                     pd[i][j][k].add_ens_var_sums(l, fcst_v);
+                  }
 
                   // Apply observation error perturbation, if requested
                   if(obs_error_info->flag) {
@@ -1781,6 +1793,21 @@ void VxPairDataEnsemble::calc_obs_summary() {
 
 ////////////////////////////////////////////////////////////////////////
 
+void VxPairDataEnsemble::set_ctrl_index(int index) {
+
+   for(int i=0; i < n_msg_typ; i++){
+      for(int j=0; j < n_mask; j++){
+         for(int k=0; k < n_interp; k++){
+            pd[i][j][k].ctrl_index = index;
+         }
+      }
+   }
+
+   return;
+}
+
+////////////////////////////////////////////////////////////////////////
+
 void VxPairDataEnsemble::set_skip_const(bool tf) {
 
    for(int i=0; i < n_msg_typ; i++){
diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.h b/met/src/libcode/vx_statistics/pair_data_ensemble.h
index e9c2417e49..bee6b11838 100644
--- a/met/src/libcode/vx_statistics/pair_data_ensemble.h
+++ b/met/src/libcode/vx_statistics/pair_data_ensemble.h
@@ -90,6 +90,7 @@ class PairDataEnsemble : public PairBase {
 
       int        n_ens;           // Number of ensemble members
       int        n_pair;          // Number of valid pairs, n_obs - sum(skip_ba)
+      int        ctrl_index;      // Index of the control member
       bool       skip_const;      // Skip cases where the observation and
                                   // all ensemble members are constant
       BoolArray  skip_ba;         // Flag for each observation [n_obs]
@@ -106,6 +107,7 @@ class PairDataEnsemble : public PairBase {
 
       NumArray   esum_na;         // Sum of unperturbed ensemble values [n_obs]
       NumArray   esumsq_na;       // Sum of unperturbed ensemble squared values [n_obs]
+      NumArray   esumn_na;        // Count of ensemble values [n_obs]
 
       NumArray   mn_na;           // Ensemble mean value [n_obs]
       NumArray   mn_oerr_na;      // Mean of perturbed members [n_obs]
@@ -293,6 +295,8 @@ class VxPairDataEnsemble {
 
       void calc_obs_summary();
 
+      void set_ctrl_index(int);
+
       void set_skip_const(bool);
 };
 

From fa5bed5985b6868533c36f65930bcbc9ee31ff49 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 16:37:04 -0600
Subject: [PATCH 13/19] Per #1905, store the index of the ensemble control
 member.

---
 met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc | 1 +
 met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h  | 3 +++
 2 files changed, 4 insertions(+)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
index 86d98605c3..da1f28c7a8 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
@@ -936,6 +936,7 @@ void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info) {
    vx_pd.set_duplicate_flag(duplicate_flag);
    vx_pd.set_obs_summary(obs_summary);
    vx_pd.set_obs_perc_value(obs_perc);
+   vx_pd.set_ctrl_index(ctrl_index);
    vx_pd.set_skip_const(skip_const);
 
    return;
diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
index 289ac667d6..eb0f18701b 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
@@ -26,6 +26,9 @@
 
 ////////////////////////////////////////////////////////////////////////
 
+// Index for the control memeber in the ensemble file list
+static int ctrl_index = bad_data_int;
+
 // Indices for the output flag types in the configuration file
 static const int i_ecnt     = 0;
 static const int i_rps      = 1;

From 2b7952a7261b630fe6d01c979d0e2aebaa4770d3 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 16:38:20 -0600
Subject: [PATCH 14/19] Per #1905, store the ensemble control as the first
 member.

---
 .../tools/core/ensemble_stat/ensemble_stat.cc | 97 +++++++++----------
 1 file changed, 48 insertions(+), 49 deletions(-)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
index fc39e16776..5ff66f6d1a 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -293,6 +293,13 @@ void process_command_line(int argc, char **argv) {
       exit(1);
    }
 
+   // Prepend the control member, if specified
+   if(ctrl_file.nonempty()) {
+      ctrl_index = 0;
+      ens_file_list.insert(ctrl_index, ctrl_file.c_str());
+      n_ens++;
+   }
+
    // Check that the end_ut >= beg_ut
    if(obs_valid_beg_ut != (unixtime) 0 &&
       obs_valid_end_ut != (unixtime) 0 &&
@@ -383,13 +390,10 @@ void process_command_line(int argc, char **argv) {
    mlog << Debug(1) << "Ensemble Files["
         << ens_file_list.n() << "]:\n";
    for(i=0; i<ens_file_list.n(); i++) {
-      mlog << "   " << ens_file_list[i]  << "\n";
+      mlog << "   " << ens_file_list[i]
+           << (i == ctrl_index ? " (control)\n" : "\n");
    }
 
-   // List the control member file
-   mlog << Debug(1) << "Ensemble Control: "
-        << (ctrl_file.empty() ? "None" : ctrl_file.c_str()) << "\n";
-
    // List the input gridded observations files
    if(grid_obs_file_list.n() > 0) {
       mlog << Debug(1) << "Gridded Observation Files["
@@ -761,8 +765,8 @@ bool get_data_plane_array(const char *infile, GrdFileType ftype,
 
 void process_ensemble() {
    int i, j;
-   bool need_reset;
-   DataPlane ens_dp, ctrl_dp;
+   bool reset;
+   DataPlane ens_dp;
    unixtime max_init_ut = bad_data_ll;
 
    // Loop through each of the ensemble fields to be processed
@@ -772,18 +776,38 @@ void process_ensemble() {
            << "Processing ensemble field: "
            << conf_info.ens_info[i]->magic_str() << "\n";
 
-      // Need to reinitialize counts and sums for each ensemble field
-      need_reset = true;
-
       // Loop through each of the input forecast files
-      for(j=0; j<ens_file_list.n(); j++) {
+      for(j=0, reset=true; j<ens_file_list.n(); j++) {
 
          // Skip bad data files
-         if(!ens_file_vld[j]) continue;
+         if(!ens_file_vld[j]) {
+
+            // Error out if the control member file is not valid
+            if(j==ctrl_index) {
+               mlog << Error << "\nprocess_ensemble() -> "
+                    << "the control member file is not valid.\n\n";
+               exit(1);
+            }
+
+            // Otherwise, continue
+            continue;
+         }
 
          // Read the current field
          if(!get_data_plane(ens_file_list[j].c_str(), etype,
-                            conf_info.ens_info[i], ens_dp, true)) continue;
+                            conf_info.ens_info[i], ens_dp, true)) {
+
+            // Error out if the control member data cannot be read
+            if(j==ctrl_index) {
+               mlog << Error << "\nprocess_ensemble() -> "
+                    << "can't find \"" << conf_info.ens_info[i]->magic_str()
+                    << "\" data in the control member file.\n\n";
+               exit(1);
+            }
+
+            // Otherwise, continue
+            continue;
+         }
 
          // Create a NetCDF file to store the ensemble output
          if(nc_out == (NcFile *) 0) {
@@ -791,40 +815,13 @@ void process_ensemble() {
          }
 
          // Reset the running sums and counts
-         if(need_reset) {
-
-            need_reset = false;
-
-            // Reset the running sums and counts
+         if(reset) {
             clear_counts();
-
-            // Read ensemble control member data, if provided
-            if(ctrl_file.nonempty()) {
-
-               // Error out if missing
-               if(!get_data_plane(ctrl_file.c_str(), etype,
-                                  conf_info.ens_info[i], ctrl_dp, true)) {
-                  mlog << Error << "\nprocess_ensemble() -> "
-                       << "control member ensemble field \""
-                       << conf_info.ens_info[i]->magic_str()
-                       << "\" not found in file \"" << ctrl_file << "\"\n\n";
-                  exit(1);
-               }
-
-               // Apply current data to the running sums and counts
-               track_counts(i, ctrl_dp, true);
-
-            } // end if ctrl_file
-
-            mlog << Debug(3)
-                 << "Found " << (ctrl_dp.is_empty() ? 0 : 1)
-                 << " control member field(s) for \""
-                 << conf_info.ens_info[i]->magic_str() << "\".\n";
-
-         } // end if need_reset
+            reset = false;
+         }
 
          // Apply current data to the running sums and counts
-         track_counts(i, ens_dp, false);
+         track_counts(i, ens_dp, j == ctrl_index);
 
          // Keep track of the maximum initialization time
          if(is_bad_data(max_init_ut) || ens_dp.init() > max_init_ut) {
@@ -838,8 +835,7 @@ void process_ensemble() {
       write_ens_nc(i, ens_dp);
 
       // Store the ensemble mean output file
-      ens_mean_file =
-         out_nc_file_list[out_nc_file_list.n() - 1];
+      ens_mean_file = out_nc_file_list[out_nc_file_list.n() - 1];
 
    } // end for i
 
@@ -861,7 +857,7 @@ void process_vx() {
       if(point_obs_file_list.n() == 0 &&
          grid_obs_file_list.n()  == 0) {
          mlog << Error << "\nprocess_vx() -> "
-              << " when \"fcst.field\" is non-empty, you must use "
+              << "when \"fcst.field\" is non-empty, you must use "
               << "\"-point_obs\" and/or \"-grid_obs\" to specify the "
               << "verifying observations.\n\n";
          exit(1);
@@ -1116,7 +1112,8 @@ int process_point_ens(int i_ens, int &n_miss) {
                                 ens_mean_file : ens_mean_user);
 
    mlog << Debug(2) << "\n" << sep_str << "\n\n"
-        << "Processing " << file_type << " file: " << ens_file << "\n";
+        << "Processing " << file_type << " file: " << ens_file
+        << (i_ens == ctrl_index ? " (control)\n" : "\n");
 
    // Loop through each of the fields to be verified and extract
    // the forecast fields for verification
@@ -1687,6 +1684,7 @@ void process_grid_vx() {
             pd_all.clear();
             pd_all.set_ens_size(n_vx_vld[i]);
             pd_all.set_climo_cdf_info(conf_info.vx_opt[i].cdf_info);
+            pd_all.ctrl_index = conf_info.vx_opt[i].vx_pd.pd[0][0][0].ctrl_index;
             pd_all.skip_const = conf_info.vx_opt[i].vx_pd.pd[0][0][0].skip_const;
 
             // Apply the current mask to the fields and compute the pairs
@@ -1916,7 +1914,8 @@ void process_grid_scores(int i_vx,
             pd.add_ens(j-n_miss, fcst_dp[j](x, y));
 
             // Store the unperturbed ensemble value
-            pd.add_ens_var_sums(i, fraw_dp[j](x, y));
+            // Exclude the control member from the variance
+            if(j != ctrl_index) pd.add_ens_var_sums(i, fraw_dp[j](x, y));
          }
       } // end for j
    } // end for i

From 2b02ea65e5694fa28bb1f3283ca549825b615183 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 17:37:04 -0600
Subject: [PATCH 15/19] Per #1905, clean up comments.

---
 met/src/libcode/vx_statistics/pair_data_ensemble.cc | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.cc b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
index 77b66536b3..0d2fc23e29 100644
--- a/met/src/libcode/vx_statistics/pair_data_ensemble.cc
+++ b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
@@ -433,7 +433,7 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
             // Compute perturbed ensemble mean and variance
             // Exclude the control member from the variance
             mn_oerr_na.add(cur_ens.mean());
-            var_oerr_na.add(cur_ens.variance(ctrl_index)); // JHG
+            var_oerr_na.add(cur_ens.variance(ctrl_index));
 
             // Compute the variance plus observation error variance.
             var_plus_oerr_na.add(var_unperturbed +
@@ -476,8 +476,8 @@ void PairDataEnsemble::compute_pair_vals(const gsl_rng *rng_ptr) {
 
          // Ensemble mean and standard deviation
          // Exclude the control member from the standard deviation
-         mean = cur_ens.mean();
-         stdev = cur_ens.stdev(ctrl_index); // JHG
+         mean  = cur_ens.mean();
+         stdev = cur_ens.stdev(ctrl_index);
 
          // Store Gaussian CRPS stats
          crps_gaus_na.add(compute_crps_gaus(o_na[i], mean, stdev));

From 872c50fe7c8585d277619d6f44e06a4592b41f5f Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 17:38:10 -0600
Subject: [PATCH 16/19] Per #1905, define ctrl_index in ensemble_stat.h instead
 of ensemble_stat_conf_info.h which was confusing the scope in which it was
 defined.

---
 met/src/tools/core/ensemble_stat/ensemble_stat.cc          | 2 +-
 met/src/tools/core/ensemble_stat/ensemble_stat.h           | 1 +
 .../tools/core/ensemble_stat/ensemble_stat_conf_info.cc    | 6 +++---
 met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h | 7 ++-----
 4 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
index 5ff66f6d1a..e36cb92a3d 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -867,7 +867,7 @@ void process_vx() {
       conf_info.process_masks(grid);
 
       // Setup the PairDataEnsemble objects
-      conf_info.set_vx_pd(n_vx_vld);
+      conf_info.set_vx_pd(n_vx_vld, ctrl_index);
 
       // Process the point observations
       if(point_obs_flag) process_point_vx();
diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.h b/met/src/tools/core/ensemble_stat/ensemble_stat.h
index 784eecad0e..c40771d1da 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.h
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.h
@@ -106,6 +106,7 @@ static ConcatString ens_mean_user; // User-specified ensemble mean data file
 static ConcatString ens_mean_file; // Computed ensemble mean output file
 
 static ConcatString ctrl_file;     // Control member
+static int          ctrl_index = bad_data_int; // Control member index
 
 // Input Observation files
 static StringArray  grid_obs_file_list;
diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
index da1f28c7a8..b338fe4f3f 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.cc
@@ -590,12 +590,12 @@ int EnsembleStatConfInfo::n_stat_row() const {
 
 ////////////////////////////////////////////////////////////////////////
 
-void EnsembleStatConfInfo::set_vx_pd(const IntArray &ens_size) {
+void EnsembleStatConfInfo::set_vx_pd(const IntArray &ens_size, int ctrl_index) {
 
    // This should be called after process_masks()
    for(int i=0; i<n_vx; i++) {
 
-      vx_opt[i].set_vx_pd(this);
+      vx_opt[i].set_vx_pd(this, ctrl_index);
 
       // Set up the ensemble size
       vx_opt[i].vx_pd.set_ens_size(ens_size[i]);
@@ -840,7 +840,7 @@ void EnsembleStatVxOpt::process_config(GrdFileType ftype, Dictionary &fdict,
 
 ////////////////////////////////////////////////////////////////////////
 
-void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info) {
+void EnsembleStatVxOpt::set_vx_pd(EnsembleStatConfInfo *conf_info, int ctrl_index) {
    int i, n;
    int n_msg_typ = msg_typ.n_elements();
    int n_mask    = mask_name.n_elements();
diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
index eb0f18701b..7148109bba 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat_conf_info.h
@@ -26,9 +26,6 @@
 
 ////////////////////////////////////////////////////////////////////////
 
-// Index for the control memeber in the ensemble file list
-static int ctrl_index = bad_data_int;
-
 // Indices for the output flag types in the configuration file
 static const int i_ecnt     = 0;
 static const int i_rps      = 1;
@@ -142,7 +139,7 @@ class EnsembleStatVxOpt {
       void process_config(GrdFileType, Dictionary &,
                           GrdFileType, Dictionary &,
                           gsl_rng *, bool, bool, bool);
-      void set_vx_pd(EnsembleStatConfInfo *);
+      void set_vx_pd(EnsembleStatConfInfo *, int);
 
       void set_perc_thresh(const PairDataEnsemble *);
 
@@ -241,7 +238,7 @@ class EnsembleStatConfInfo {
       void process_flags ();
       void parse_nc_info ();
       void process_masks (const Grid &);
-      void set_vx_pd     (const IntArray &);
+      void set_vx_pd     (const IntArray &, int);
 
       // Dump out the counts
       int get_n_ens_var()    const;

From 07e7d516a36f79aa13ac2b43ffbe530f5156c1b2 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Fri, 29 Oct 2021 18:01:59 -0600
Subject: [PATCH 17/19] Per #1905, correct ssvar logic for skipping the
 ensemble control member.

---
 met/src/libcode/vx_statistics/pair_data_ensemble.cc | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/met/src/libcode/vx_statistics/pair_data_ensemble.cc b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
index 0d2fc23e29..4811f66147 100644
--- a/met/src/libcode/vx_statistics/pair_data_ensemble.cc
+++ b/met/src/libcode/vx_statistics/pair_data_ensemble.cc
@@ -652,11 +652,14 @@ void PairDataEnsemble::compute_ssvar() {
    for(i=0; i<o_na.n(); i++) {
 
       // Check if point should be skipped
-      // Exclude the control member from the variance
-      if(skip_ba[i] || i == ctrl_index) continue;
+      if(skip_ba[i]) continue;
 
       // Store ensemble values for the current observation
-      for(j=0, cur.erase(); j<n_ens; j++) cur.add(e_na[j][i]);
+      // Exclude the control member from the variance
+      for(j=0, cur.erase(); j<n_ens; j++) {
+         if(j == ctrl_index) continue;
+         cur.add(e_na[j][i]);
+      }
 
       // Build a variance point
       ens_ssvar_pt pt;

From 2992759d801e771a7d7975572bd875a273a69562 Mon Sep 17 00:00:00 2001
From: John Halley Gotway <johnhg@ucar.edu>
Date: Mon, 1 Nov 2021 13:12:42 -0600
Subject: [PATCH 18/19] Per #1905, revert variable name usage to minimize diffs
 from develop branch.

---
 .../tools/core/ensemble_stat/ensemble_stat.cc  | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/met/src/tools/core/ensemble_stat/ensemble_stat.cc b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
index e36cb92a3d..783691945b 100644
--- a/met/src/tools/core/ensemble_stat/ensemble_stat.cc
+++ b/met/src/tools/core/ensemble_stat/ensemble_stat.cc
@@ -2015,7 +2015,7 @@ void clear_counts() {
 
 void track_counts(int i_vx, const DataPlane &ens_dp, bool is_ctrl) {
    int i, j, k;
-   double ens;
+   double v;
 
    // Ensemble thresholds
    const int n_thr = conf_info.ens_ta[i_vx].n();
@@ -2025,10 +2025,10 @@ void track_counts(int i_vx, const DataPlane &ens_dp, bool is_ctrl) {
    for(i=0; i<nxy; i++) {
 
       // Get current values
-      ens = ens_dp.data()[i];
+      v = ens_dp.data()[i];
 
       // Skip bad data values
-      if(is_bad_data(ens)) continue;
+      if(is_bad_data(v)) continue;
 
       // Otherwise, update counts and sums
       else {
@@ -2037,22 +2037,22 @@ void track_counts(int i_vx, const DataPlane &ens_dp, bool is_ctrl) {
          cnt_na.buf()[i] += 1;
 
          // Ensemble sum
-         sum_na.buf()[i] += ens;
+         sum_na.buf()[i] += v;
 
          // Ensemble min and max
-         if(ens <= min_na.buf()[i] || is_bad_data(min_na.buf()[i])) min_na.buf()[i] = ens;
-         if(ens >= max_na.buf()[i] || is_bad_data(max_na.buf()[i])) max_na.buf()[i] = ens;
+         if(v <= min_na.buf()[i] || is_bad_data(min_na.buf()[i])) min_na.buf()[i] = v;
+         if(v >= max_na.buf()[i] || is_bad_data(max_na.buf()[i])) max_na.buf()[i] = v;
 
          // Standard deviation sum, sum of squares, and count, excluding control member
          if(!is_ctrl) {
-            stdev_sum_na.buf()[i] += ens;
-            stdev_ssq_na.buf()[i] += ens*ens;
+            stdev_sum_na.buf()[i] += v;
+            stdev_ssq_na.buf()[i] += v*v;
             stdev_cnt_na.buf()[i] += 1;
          }
 
          // Event frequency
          for(j=0; j<n_thr; j++) {
-            if(thr_buf[j].check(ens)) thresh_cnt_na[j].inc(i, 1);
+            if(thr_buf[j].check(v)) thresh_cnt_na[j].inc(i, 1);
          }
       } // end else
    } // end for i

From a616008acd0ccb67e048649328e437f0d8769c7d Mon Sep 17 00:00:00 2001
From: j-opatz <59586397+j-opatz@users.noreply.github.com>
Date: Tue, 2 Nov 2021 13:30:59 -0600
Subject: [PATCH 19/19] Small spelling correction

---
 met/docs/Users_Guide/ensemble-stat.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/met/docs/Users_Guide/ensemble-stat.rst b/met/docs/Users_Guide/ensemble-stat.rst
index cf2cbd0012..1c131a9580 100644
--- a/met/docs/Users_Guide/ensemble-stat.rst
+++ b/met/docs/Users_Guide/ensemble-stat.rst
@@ -97,7 +97,7 @@ Optional arguments for ensemble_stat
 
 6. To override the simple ensemble mean value of the input ensemble members for the ECNT, SSVAR, and ORANK line types, the **-ens_mean file** option specifies an ensemble mean model data file. This option replaces the **-ssvar_mean file** option from earlier versions of MET.
 
-7. The **-ctrl file** option specifies and ensemble control member data file. The control member is included in the computation of the ensemble mean but excluded from the spread.
+7. The **-ctrl file** option specifies an ensemble control member data file. The control member is included in the computation of the ensemble mean but excluded from the spread.
 
 8. To filter point observations by time, use **-obs_valid_beg time** in YYYYMMDD[_HH[MMSS]] format to set the beginning of the matching observation time window.