// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
// ** Copyright UCAR (c) 1992 - 2007
// ** University Corporation for Atmospheric Research (UCAR)
// ** National Center for Atmospheric Research (NCAR)
// ** Research Applications Lab (RAL)
// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*

///////////////////////////////////////////////////////////////////////
//
//   Filename:   grid_stat.cc
//
//   Description:
//
//   Mod#   Date      Name            Description
//   ----   ----      ----            -----------
//   000    04/18/07  Halley Gotway   New
//   001    10/05/07  Halley Gotway   Add checks to only compute
//          confidence intervals for sufficiently large sample sizes.
//   002    10/19/07  Halley Gotway   Add ability to smooth the forecast
//          field using the interp_method and interp_width parameters
//   003    02/12/08  Halley Gotway   Fix bug in writing COP line to
//                    write out OY and ON rather than FY and FN.
//
////////////////////////////////////////////////////////////////////////

using namespace std;

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <ctype.h>
#include <dirent.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>

#include "netcdf.hh"
#include "vx_grib_classes/grib_classes.h"

#include "grid_stat_Conf.h"

#include "vx_wrfdata/vx_wrfdata.h"
#include "vx_met_util/vx_met_util.h"
#include "vx_data_grids/grid.h"
#include "vx_util/vx_util.h"
#include "vx_cal/vx_cal.h"
#include "vx_math/vx_math.h"
#include "vx_contable/vx_contable.h"
#include "vx_gsl_prob/vx_gsl_prob.h"
#include "vx_econfig/result.h"

////////////////////////////////////////////////////////////////////////
//
// Create output files using the following naming convention:
//
//    grid_stat_HHMMSSL_YYYYMMDD_HHMMSSV.vsdb/.txt/.nc
//
//    Where L indicates lead time and V indicates valid time.
//
////////////////////////////////////////////////////////////////////////

static const char * program_name = "grid_stat";

// Input configuration file
static char config_file[max_str_len];

// Input forecast Grib or NetCDF files
static char fcst_file[max_str_len];
static GribFile fcst_gb_file;
static NcFile  *fcst_nc_file; 

// Input observation Grib or NetCDF files
static char obs_file[max_str_len];
static GribFile obs_gb_file;
static NcFile  *obs_nc_file; 

// Output NetCDF file
static char    out_nc_file[max_str_len];
static NcFile  *nc_out = (NcFile *) 0;
static NcDim   *lat_dim, *lon_dim;
static NcVar   *lat_var, *lon_var;

static const int max_num_vars = 50;
static const int var_str_len  = 32;
static const int max_line_len = 1024;
// Exponent used in distance weighted mean calculations
static const int i_pow = 2;

// Grid definition data structures
static PlateCarreeData   pc_data;
static LambertData       lc_data;
static StereographicData st_data;

// Variables for command line arguments
static grid_stat_Conf conf;
static char out_dir[max_str_len];
static int verbosity = 1;

// Indices for the output flag types in the configuration file
static const int i_fho   = 0;
static const int i_ctc   = 1;
static const int i_ctp   = 2;
static const int i_cfp   = 3;
static const int i_cop   = 4;
static const int i_cts   = 5;
static const int i_cnt   = 6;
static const int i_sl1l2 = 7;
static const int i_nc    = 8;

// Grid variables
static Grid grid, fcst_grid, obs_grid;

// Flags to indicate the format of the input files. 
//    0 = netCDF (output of pcp_combine)
//    1 = Grib Version 1
static int fcst_fmt_flag = 0;
static int obs_fmt_flag = 0;

// Strings to be output in the VSDB file
static const char *met_version = "V01";
static char model_name[max_str_len];
static char vsdb_code[max_str_len];
static char field_name[max_str_len];
static char level_name[max_str_len];

// Number of Grib Codes to be verified
static int n_gc;

// Objects to hold Grib Code information
// Sized as gc_info[n_gc]
static GCInfo *gc_info = (GCInfo *) 0;

// Objects to hold threshold information for each Grib Code
// Sized as gc_ti[n_gc]
static ThreshInfo *gc_ti = (ThreshInfo *) 0;

// Multiple smoothing operations to be performed on the forecast field
// prior to verification
static int  n_interp;
static int *interp_mthd = (int *) 0;
static int *interp_wdth = (int *) 0;

// Multiple verification masking regions to be applied to the forecast
// and observation fields
static int      n_masks;
static WrfData *mask_wd    = (WrfData *) 0;
static char   **mask_names = (char **) 0;

// Multiple values of alpha to be used when computing confidence intervals
// for continuous variables
static int n_ci_alpha;

// Class to store Contingency Table Counts and statistics
static int n_cts;
static CTSInfo **cts_info = (CTSInfo **) 0;

// Class to store Continuous Statistics
static int n_cnt;
static CNTInfo **cnt_info = (CNTInfo **) 0;

// Output file streams for creating VSDB and text files
// 
static ofstream *vsdb_out, *cts_out, *cnt_out, *sl1l2_out;
static ofstream *fho_out, *ctc_out, *ctp_out, *cfp_out, *cop_out;
static char vsdb_file[max_str_len];
static char fho_file[max_str_len], ctc_file[max_str_len];
static char ctp_file[max_str_len], cfp_file[max_str_len], cop_file[max_str_len];
static char cts_file[max_str_len], cnt_file[max_str_len], sl1l2_file[max_str_len]; 

///////////////////////////////////////////////////////////////////////////////

static void process_command_line(int, char **);
static void process_fcst_obs_files();
static void process_config();
static void process_fields();
static void process_masks();
static void clean_up();

static void setup_first_pass(const WrfData &);
static void setup_txt_files(unixtime, int);
static void setup_nc_file(unixtime, int);

static void compute_cts(const WrfData &, const WrfData &, int, int);
static void compute_cnt(const WrfData &, const WrfData &, int, int);

static void write_txt_files(const WrfData &, const WrfData &, int, int);
static void write_fcst_netcdf(const WrfData &, const WrfData &, int, int, int);
static void write_obs_netcdf(const WrfData &, int);
static void close_txt_files();

static int  get_interp_mthd(const char *);
static WrfData smooth_field(WrfData *, int, int);

static void usage(int, char **);

////////////////////////////////////////////////////////////////////////

int main(int argc, char *argv[]) {

   // Process the command line arguments
   process_command_line(argc, argv);

   // Process the forecast and observation input files
   process_fcst_obs_files();

   // Process the configuration file
   process_config();   
   
   // Process the fields to be verified
   process_fields();
   
   // Close the text files and deallocate memory
   clean_up();
     
   return(0);
}   

////////////////////////////////////////////////////////////////////////

void process_command_line(int argc, char **argv) {
   int i;

   sprintf(out_dir, "%s/out/grid_stat", MET_BASE);
   
   if(argc < 4) {
      usage(argc, argv);
      exit(1);
   }

   // Store the input forecast and observation file names
   strcpy(fcst_file, argv[1]);
   strcpy(obs_file, argv[2]);
   strcpy(config_file, argv[3]);

   // Parse command line arguments
   for(i=0; i<argc; i++) {
   
      if(strcmp(argv[i], "-outdir") == 0) {
         strcpy(out_dir, argv[i+1]);
         i++;
      }
      else if(strcmp(argv[i], "-v") == 0) {
         verbosity = atoi(argv[i+1]);
         i++;
      }
   }
   
   // Read the config file
   conf.read(config_file);
   
   // Store the model name
   strcpy(model_name, conf.model().sval());

   // List the input files
   if(verbosity > 0) {
      cout << "Forecast File: " << fcst_file << "\n" 
           << "Observation File: " << obs_file << "\n" 
           << "Configuration File: " << config_file << "\n" << flush;
   }
   
   return;
}
////////////////////////////////////////////////////////////////////////

void process_fcst_obs_files() {

   // Attempt to open the forecast file as a NetCDF file.  If it
   // fails, attempt to read it as Grib.   
   fcst_nc_file = new NcFile(fcst_file);
   
   if(!fcst_nc_file || !(fcst_nc_file->is_valid())) {
      fcst_nc_file->close();
      delete fcst_nc_file;
      fcst_nc_file = (NcFile *) 0;
      fcst_fmt_flag = 1;
   }
   else {
      fcst_fmt_flag = 0;
   }

   // Attempt to open the forecast file as a Grib file.
   if(fcst_fmt_flag == 1) {

      if( !(fcst_gb_file.open(fcst_file)) ) {
         cerr << "\n\nERROR: process_fcst_obs_files() -> "
              << "can't open grib file: "
              << fcst_file << "\n\n" << flush;
         exit(1);
      }
   }
   
   // Attempt to open the observation file as a NetCDF file.  If it
   // fails, attempt to read it as Grib.   
   obs_nc_file = new NcFile(obs_file);
   
   if(!obs_nc_file || !(obs_nc_file->is_valid())) {
      obs_nc_file->close();
      delete obs_nc_file;
      obs_nc_file = (NcFile *) 0;
      obs_fmt_flag = 1;
   }
   else {
      obs_fmt_flag = 0;
   }

   // Attempt to open the observation file as a Grib file.
   if(obs_fmt_flag == 1) {

      if(!(obs_gb_file.open(obs_file)) ) {
         cerr << "\n\nERROR: process_fcst_obs_files() -> "
              << "can't open grib file: "
              << obs_file << "\n\n" << flush;
         exit(1);
      }
   }
   
   return;
}

////////////////////////////////////////////////////////////////////////

void process_config() {
   int i, j, n, wdth, n_mthd, n_wdth;
   
   //
   // Conf: vx_grib_code
   //
   
   // Parse out the grib codes to be verified
   n_gc = conf.n_vx_grib_code_elements();
   
   if(n_gc == 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one value must be provided "
           << "for vx_grib_code.\n\n" << flush;      
      exit(1);
   }
      
   // Allocate space to store the grib code information
   gc_info = new GCInfo [n_gc];
      
   for(i=0; i<n_gc; i++) {
      parse_gcinfo(conf.vx_grib_code(i).sval(), gc_info[i]);
      
      // No support for WDIR
      if(gc_info[i].code == wdir_grib_code) {
         cerr << "\n\nERROR: process_config() -> "
              << "the wind direction field may not be verified using grid_stat.\n\n"
              << flush;
         exit(1);
      }
   }

   // Grib File Input
   // If fcst_fmt_flag == 1 or obs_fmt_flag == 1, check that at least
   // one grib code is supplied
   if( (fcst_fmt_flag == 1 || obs_fmt_flag == 1) && n_gc <= 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one gc_info must be provided "
           << "when processing grib input files.\n\n" << flush;
      exit(1);
   }   
      
   // NetCDF File Input
   // If mixing NetCDF with Grib input, check that the grib code specified
   // is accumulated precip (APCP)
   if( (fcst_fmt_flag == 0 || obs_fmt_flag == 0)
    && (n_gc != 1 || gc_info[0].code != apcp_grib_code) ) {
      cerr << "\n\nERROR: process_config() -> " 
           << "When using NetCDF input, gc_info must "
           << "contain only grib code, " << apcp_grib_code << ".\n\n" 
           << flush;
      exit(1);
   }
   
   //
   // Conf: thresholds
   //
   
   // Check that the number of forecast threshold levels matches n_gc
   if(conf.n_thresholds_elements() != n_gc) {
      cerr << "\n\nERROR: process_config() -> " 
           << "The number threshold levels provided must match the number "
           << "of grib codes provided in vx_grib_code.\n\n"
           << flush;
      exit(1);
   }
   
   // Allocate space to store the threshold information
   gc_ti = new ThreshInfo [n_gc];
   
   // Parse the threshold information
   for(i=0; i<n_gc; i++)
      gc_ti[i].parse_thresh_str(conf.thresholds(i).sval());

   //
   // Conf: ci_alpha
   //
   
   // Check that at least one alpha value is provided
   if((n_ci_alpha = conf.n_ci_alpha_elements()) == 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one confidence interval alpha value must be "
           << "specified.\n\n" << flush;
      exit(1);
   }
   
   // Check that the values for alpha are between 0 and 1
   for(i=0; i<n_ci_alpha; i++) {
      if(conf.ci_alpha(i).dval() <= 0.0 || conf.ci_alpha(i).dval() >= 1.0) {
         cerr << "\n\nERROR: process_config() -> " 
              << "All confidence interval alpha values (" 
              << conf.ci_alpha(i).dval() << ") must be greater than 0 "
              << "and less than 1.\n\n" << flush;
         exit(1);
      }
   }
   
   //
   // Conf: interp_method
   //
   
   // Check that at least one interpolation method is provided
   if((n_mthd = conf.n_interp_method_elements()) <= 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one interpolation method must be provided.\n\n"
           << flush;
      exit(1);
   }

   // Check that distance-weighted mean and least squares fit are not selected
   for(i=0; i<n_mthd; i++) {
      if(get_interp_mthd(conf.interp_method(i).sval()) == i_dw_mean ||
         get_interp_mthd(conf.interp_method(i).sval()) == i_ls_fit) {
         cerr << "\n\nERROR: process_config() -> " 
              << "The interpolation method may not be set to DW_MEAN "
              << "or LS_FIT for grid_stat.\n\n" << flush;
         exit(1);
      }
   }
   
   //
   // Conf: interp_width
   //   

   // Check that at least one interpolation width is provided
   if((n_wdth = conf.n_interp_width_elements()) <= 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one interpolation width must be provided.\n\n"
           << flush;
      exit(1);
   }
      
   // Do error checking and compute the total number of 
   // interpolations to be performed
   n_interp = 0;
   for(i=0; i<n_wdth; i++) {
   
      if(conf.interp_width(i).ival() < 1 ||
         conf.interp_width(i).ival()%2 == 0) {
         cerr << "\n\nERROR: process_config() -> " 
              << "The interpolation width must be set "
              << "to odd values greater than or equal to 1 (" 
              << conf.interp_width(i).ival() << ").\n\n" << flush;
         exit(1);
      }
      
      wdth = conf.interp_width(i).ival();
      if(wdth == 1) n_interp += 1;
      if(wdth >  1) n_interp += n_mthd;
   }
   
   // Allocate space for the interpolation methods and widths
   interp_mthd = new int [n_interp];
   interp_wdth = new int [n_interp];
   
   // Set each interpolation method and width
   n = 0;
   for(i=0; i<n_wdth; i++) {
   
      interp_wdth[n] = conf.interp_width(i).ival();
      
      // For an interpolation width of 1, set the method the method to
      // unweighted mean - which is really just nearest neighbor
      if(interp_wdth[n] == 1) {
         interp_mthd[n] = i_uw_mean;
         n++;
         continue;
      }
   
      for(j=0; j<n_mthd; j++) {
         interp_mthd[n] = get_interp_mthd(conf.interp_method(j).sval());
         n++;
      }
   }

   //
   // Conf: interp_threshold
   //
   
   // Check that the interpolation threshold is set between 
   // 0 and 1.
   if(conf.interp_threshold().dval() < 0.0 
   || conf.interp_threshold().dval() > 1.0) {
         cerr << "\n\nERROR: process_config() -> " 
              << "The interpolation threshold value must be set "
              << "between 0 and 1.\n\n" << flush;
         exit(1);
   }   
      
   //
   // Conf: output_flag
   //
   
   // Check that at least one output VSDB type is requested
   for(i=i_fho, n=0; i<i_sl1l2; i++) n += (conf.output_flag(i).ival() > 0.0);
   if(n == 0) {
      cerr << "\n\nERROR: process_config() -> " 
           << "At least one output VSDB type must be requested.\n\n" 
           << flush;
      exit(1);
   }

   return;
}

////////////////////////////////////////////////////////////////////////

void process_fields() {
   int i, j, status;
   WrfData fcst_wd, obs_wd, fcst_wd_smooth;
   GribRecord fcst_r, obs_r;

   // Loop through each of the fields to be verified
   for(i=0; i<n_gc; i++) {

      if(verbosity > 1 && gc_info[i].lvl_str != (char *) 0) {
         cout << "\n----------------------------------------\n\n"
              << "Processing Grib Code " << gc_info[i].code 
              << " with accumulation/level indicator of " << gc_info[i].lvl_str
              << "\n" << flush;
      }
      else if(verbosity > 1) {
         cout << "\n----------------------------------------\n\n"
              << "Processing Grib Code " << gc_info[i].code
              << " with no accumulation/level indicator\n" << flush;
      }
   
      // Process the forecast file in NetCDF format (output of pcp_combine)
      if(fcst_fmt_flag == 0) {
         read_pcp_combine_netcdf(fcst_nc_file, fcst_wd, fcst_grid, verbosity, 
                                 pc_data, lc_data, st_data);
      }
      // Process the forecast file in Grib format
      else {
      
         status = get_grib_record(fcst_gb_file, fcst_r, 
                                  gc_info[i].code, gc_info[i].lvl_1, gc_info[i].lvl_2, 
                                  fcst_wd, fcst_grid, verbosity, 
                                  pc_data, lc_data, st_data);

         if(status != 0) {
            cout << "***WARNING***: process_fields() -> "
                 << "grib code " << gc_info[i].code
                 << " with accumulation/level indicator of "
                 << gc_info[i].lvl_str << " not found in grib file: "
                 << fcst_file << "\n" << flush;
            continue;
         }
      }

      // Process the observation file in NetCDF format (output of pcp_combine)
      if(obs_fmt_flag == 0) {
         read_pcp_combine_netcdf(obs_nc_file, obs_wd, obs_grid, verbosity, 
                                 pc_data, lc_data, st_data);
      }
      // Process the observation file in Grib format
      else {
      
         status = get_grib_record(obs_gb_file, obs_r, 
                                  gc_info[i].code, gc_info[i].lvl_1, gc_info[i].lvl_2, 
                                  obs_wd, obs_grid, verbosity, 
                                  pc_data, lc_data, st_data);

         if(status != 0) {
            cout << "***WARNING***: process_fields() -> "
                 << "grib code " << gc_info[i].code
                 << " with accumulation/level indicator of "
                 << gc_info[i].lvl_str << " not found in grib file: "
                 << obs_file << "\n" << flush;
            continue;
         }
      }

      // Check that the grid dimensions match
      if(fcst_grid.nx() != obs_grid.nx() ||
         fcst_grid.ny() != obs_grid.ny()) {
         cerr << "\n\nERROR: process_fields() -> "
              << "Forecast and observation grid dimensions "
              << "do not match (" << fcst_grid.nx() << ", "
              << fcst_grid.ny() << ") != (" << obs_grid.nx()
              << ", " << obs_grid.ny() << ") for grib code "
              << gc_info[i].code << ".\n\n" << flush;
         exit(1);
      }
      
      // Check that the valid times match
      if(fcst_wd.get_valid_time() != obs_wd.get_valid_time()) {
         cout << "***WARNING***: process_fields() -> "
              << "Forecast and observation valid times do not match "
              << fcst_wd.get_valid_time() << " != "
              << obs_wd.get_valid_time() << " for grib code "
              << gc_info[i].code << ".\n\n" << flush;
      }
      
      // Check that the accumulation intervals match
      if(fcst_wd.get_accum_time() != obs_wd.get_accum_time()) {
         cout << "***WARNING***: process_fields() -> "
              << "Forecast and observation accumulation times do not match "
              << fcst_wd.get_accum_time() << " != "
              << obs_wd.get_accum_time() << " for grib code "
              << gc_info[i].code << ".\n\n" << flush;
      }
      
      // This is the first pass through the loop and grid is unset
      if(grid.nx() == 0 && grid.ny() == 0) {
      
         // Setup the first pass through the data
         setup_first_pass(fcst_wd);
      }   
      // For multiple verification fields, check to make sure that the
      // grid dimensions don't change
      else {
      
         if(fcst_grid.nx() != grid.nx() ||
            fcst_grid.ny() != grid.ny() ||
            obs_grid.nx()  != grid.nx() ||
            obs_grid.ny()  != grid.ny()) {
            cerr << "\n\nERROR: process_fields() -> "
                 << "The grid dimensions must remain constant ("
                 << grid.nx() << ", " << grid.ny() << ") != (" 
                 << fcst_grid.nx() << ", " << fcst_grid.ny() << ") or (" 
                 << obs_grid.nx() << ", " << obs_grid.ny() << ")\n\n" 
                 << flush;
            exit(1);
         }  
      }      

      // Setup strings for accumulated precip (APCP)
      if(fcst_fmt_flag == 0 || 
         obs_fmt_flag == 0 ||
         gc_info[i].code == apcp_grib_code) {

         strcpy(vsdb_code, "MC_PCP");
         sprintf(field_name, "APCP/%.2i", nint(fcst_wd.get_accum_time()/sec_per_hour));
         strcpy(level_name, "SFC");
      }
      // Setup strings for other grib codes
      else{
         strcpy(vsdb_code, "ANALYS");
         if(nint(fcst_wd.get_accum_time()/sec_per_hour) > 0) {
            sprintf(field_name, "%s/%.2i", 
                    get_grib_code_abbr(gc_info[i].code, 
                                       conf.ncep_defaults().ival(), status),
                    nint(fcst_wd.get_accum_time()/sec_per_hour));
         }
         else {
            strcpy(field_name, 
                   get_grib_code_abbr(gc_info[i].code, 
                                      conf.ncep_defaults().ival(), status));
         }  
         strcpy(level_name, 
                get_grib_level_str(fcst_r.pds->type, conf.ncep_defaults().ival(), 
                                   fcst_r.pds->level_info, status));
      }
      
      // Loop through and apply each of the smoothing operations
      for(j=0; j<n_interp; j++) {
      
         // Smooth the forecast field based on the interp_mthd and interp_wdth
         fcst_wd_smooth = smooth_field(&fcst_wd, interp_mthd[j], interp_wdth[j]);
      
         // Compute contingency table counts for each 
         // grib code/smoothing operation/threshold/masking region
         // if FHO, CTC, CTP, CFP, COP, or CTS is requested in the config file
         if(conf.output_flag(i_fho).ival() || conf.output_flag(i_ctc).ival() || 
            conf.output_flag(i_ctp).ival() || conf.output_flag(i_cfp).ival() || 
            conf.output_flag(i_cop).ival() || conf.output_flag(i_cts).ival()) 
            compute_cts(fcst_wd_smooth, obs_wd, j, i);

         // Compute continuous statistics for each 
         // grib code/smoothing operation/masking region 
         // if CNT or SL1L2 is requested in the config file
         if(conf.output_flag(i_cnt).ival() || conf.output_flag(i_sl1l2).ival()) 
            compute_cnt(fcst_wd_smooth, obs_wd, j, gc_info[i].code);

         // Create output VSDB statistics file
         write_txt_files(fcst_wd_smooth, obs_wd, j, i);
      
         // Write out the fcst and difference fields in netCDF format for 
         // each grib code if requested in the config file
         if(conf.output_flag(i_nc).ival())
            write_fcst_netcdf(fcst_wd_smooth, obs_wd, i,
                              interp_mthd[j], interp_wdth[j]);
      } // end for j
      
      // Write out the observation fields in netCDF format for each grib code
      // if requested in the config file
      if(conf.output_flag(i_nc).ival())
         write_obs_netcdf(obs_wd, i);
      
   } // end for i
   
   if(verbosity > 1) {
      cout << "\n----------------------------------------\n\n" << flush;
   }
   
   return;
}

////////////////////////////////////////////////////////////////////////

void process_masks() {
   int i, n_mask_grids, n_mask_polys;
   
   // Get the number of masks
   n_mask_grids = conf.n_mask_grids_elements();
   n_mask_polys = conf.n_mask_polys_elements();
   n_masks = n_mask_grids + n_mask_polys;
   
   // Check that at least one verification masking region is provided
   if(n_masks == 0) {
      cerr << "\n\nERROR: process_masks() -> " 
           << "At least one grid or polyline verification masking "
           << "region must be provided.\n\n" << flush;
      exit(1);
   }
   
   // Allocate space to store the masking information
   mask_wd    = new WrfData [n_masks];
   mask_names = new char *[n_masks];
   
   // Parse out the masking grids
   for(i=0; i<n_mask_grids; i++)
      parse_grid_mask(conf.mask_grids(i).sval(), grid, 
                      mask_wd[i], mask_names[i]);
   
   // Parse out the masking polys
   for(i=0; i<n_mask_polys; i++)
      parse_poly_mask(conf.mask_polys(i).sval(), grid,
                      mask_wd[i+n_mask_grids], mask_names[i+n_mask_grids]);
   
   return;
}

////////////////////////////////////////////////////////////////////////

void clean_up() {
   int i, j;

   // Close the output text files that were open for writing
   close_txt_files();
   
   // Close the output NetCDF file as long as it was opened
   if(nc_out && conf.output_flag(i_nc).ival()) {

      // List the NetCDF file after it is finished
      if(verbosity > 0) {
         cout << "Output NetCDF file: "
              << out_nc_file << "\n" << flush;
      }
      nc_out->close(); 
      delete nc_out; 
      nc_out = (NcFile *) 0;
   }
   
   // Close the forecast file
   if(fcst_fmt_flag == 0) fcst_nc_file->close();
   else                   fcst_gb_file.close();

   // Close the observation file
   if(obs_fmt_flag == 0)  obs_nc_file->close();
   else                   obs_gb_file.close();
   
   // Clean up allocated memory
   if(gc_info)          { delete [] gc_info;       gc_info = (GCInfo *) 0; }
   if(gc_ti)            { delete [] gc_ti;         gc_ti = (ThreshInfo *) 0; }
   if(interp_mthd)      { delete [] interp_mthd;   interp_mthd = (int *) 0; }
   if(interp_wdth)      { delete [] interp_wdth;   interp_wdth = (int *) 0; }
   if(mask_wd)          { delete [] mask_wd;       mask_wd = (WrfData *) 0; }
   
   for(i=0; i<n_masks; i++) {
      if(mask_names[i]) { delete [] mask_names[i]; mask_names[i] = (char *) 0; }
   }
   if(mask_names)       { delete [] mask_names;    mask_names = (char **) 0; }
   
   // Deallocate CTSInfo objects
   for(i=0; i<n_interp; i++) {
      for(j=0; j<n_cts; j++) {
         cts_info[i][j].clear();
      }
   }
   for(i=0; i<n_interp; i++) {
      if(cts_info[i])   { delete [] cts_info[i];   cts_info[i] = (CTSInfo *) 0; }
   }
   if(cts_info)         { delete [] cts_info;      cts_info = (CTSInfo **) 0; }

   // Deallocate CNTInfo objects
   for(i=0; i<n_interp; i++) {
      for(j=0; j<n_cnt; j++) {
         cnt_info[i][j].clear();
      }
   }
   for(i=0; i<n_interp; i++) {
      if(cnt_info[i])   { delete [] cnt_info[i];   cnt_info[i] = (CNTInfo *) 0; }
   }
   if(cnt_info)         { delete [] cnt_info;      cnt_info = (CNTInfo **) 0; }   
   
   return;
}

////////////////////////////////////////////////////////////////////////

void setup_first_pass(const WrfData &wd) {
   int i, j, k, max_n_thresh;
      
   // Store the grid to be used through the verification
   grid = fcst_grid;

   // Process masks Grids and Polylines in the config file
   process_masks();

   // Figure out the maximum number of thresholds being applied
   // across Grib Codes
   max_n_thresh = 0;
   for(j=0; j<n_gc; j++) {
      if(gc_ti[j].n_thresh > max_n_thresh) 
         max_n_thresh = gc_ti[j].n_thresh;
   }
   
   // Allocate space for 2X2 contingency table objects to
   // store the raw counts (max_n_thresh*n_masks).
   cts_info = new CTSInfo * [n_interp];
   n_cts = max_n_thresh*n_masks;
   for(j=0; j<n_interp; j++) cts_info[j] = new CTSInfo [n_cts];

   // Allocate space for objects to store continuous
   // statistics for grib variables
   cnt_info = new CNTInfo * [n_interp];
   n_cnt = n_masks;
   for(j=0; j< n_interp; j++) cnt_info[j] = new CNTInfo [n_cnt];
   
   // Allocate enough space in each obiect and store data for 
   // all of the alpha values specified
   for(i=0; i<n_interp; i++) {
      for(j=0; j<n_cts; j++) {
         cts_info[i][j].allocate_n_alpha(n_ci_alpha);
         
         for(k=0; k<n_ci_alpha; k++)
            cts_info[i][j].alpha[k] = conf.ci_alpha(k).dval();
      } // end for j
   } // end for i

   for(i=0; i<n_interp; i++) {
      for(j=0; j<n_cnt; j++) {
         cnt_info[i][j].allocate_n_alpha(n_ci_alpha);
         
         for(k=0; k<n_ci_alpha; k++)
            cnt_info[i][j].alpha[k] = conf.ci_alpha(k).dval();
      } // end for j
   } // end for i

   // Create output text files as requested in the config file
   setup_txt_files(wd.get_valid_time(), wd.get_lead_time());

   // If requested, create a NetCDF file to store the matched pairs and
   // difference fields for each grib code and masking region
   if(conf.output_flag(i_nc).ival()) {
      setup_nc_file(wd.get_valid_time(), wd.get_lead_time());
   }

   return;
}

////////////////////////////////////////////////////////////////////////

void setup_txt_files(unixtime valid_ut, int lead_sec) {
   int mon, day, yr, hr, min, sec;
   int l_hr, l_min, l_sec;
   char header[max_line_len];

   // Initialize output file streams
   vsdb_out = fho_out = cts_out = cnt_out = sl1l2_out = (ofstream *) 0;
   ctc_out  = ctp_out = cfp_out = cop_out = (ofstream *) 0;
   
   // Create output VSDB file names for the VSDB and optional text files
   sec_to_hms(lead_sec, l_hr, l_min, l_sec);
   unix_to_mdyhms(valid_ut, mon, day, yr, hr, min, sec);
   sprintf(vsdb_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV.vsdb",
          out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(fho_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_fho.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(ctc_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_ctc.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(ctp_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_ctp.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(cfp_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_cfp.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(cop_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_cop.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(cts_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_cts.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(cnt_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_cnt.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(sl1l2_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_sl1l2.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
         
   // Create and open the VSDB output file stream
   vsdb_out = new ofstream;
   vsdb_out->open(vsdb_file);

   if(!vsdb_out) {
      cerr << "\n\nERROR: setup_txt_files() -> "
           << "unable to open output VSDB file \""
           << vsdb_file << "\"\n\n";
      exit(1);
   }
   vsdb_out->setf(ios::fixed);

   // List VSDB file as it is being created
   if(verbosity > 1) {
      cout << "Creating Output VSDB file: "
           << vsdb_file << "\n" << flush;
   }

   // Create the FHO text file
   if(conf.output_flag(i_fho).ival() >= 2) {

      // Create and open the FHO output file stream
      fho_out = new ofstream;
      fho_out->open(fho_file);

      if(!fho_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output FHO file \""
              << fho_file << "\"\n\n";
         exit(1);
      }
      fho_out->setf(ios::fixed);

      // List FHO file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output FHO file: "
              << fho_file << "\n" << flush;
      }
      
      fho_hdr_str(header);
      *fho_out << header << "\n" << flush; 
   }
   
   // Create the CTC text file
   if(conf.output_flag(i_ctc).ival() >= 2) {

      // Create and open the CTC output file stream
      ctc_out = new ofstream;
      ctc_out->open(ctc_file);

      if(!ctc_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output CTC file \""
              << ctc_file << "\"\n\n";
         exit(1);
      }
      ctc_out->setf(ios::fixed);

      // List CTC file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output CTC file: "
              << ctc_file << "\n" << flush;
      }

      ctc_hdr_str(header);
      *ctc_out << header << "\n" << flush; 
   }

   // Create the CTP text file
   if(conf.output_flag(i_ctp).ival() >= 2) {

      // Create and open the CTP output file stream
      ctp_out = new ofstream;
      ctp_out->open(ctp_file);

      if(!ctp_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output CTP file \""
              << ctp_file << "\"\n\n";
         exit(1);
      }
      ctp_out->setf(ios::fixed);

      // List CTP file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output CTP file: "
              << ctp_file << "\n" << flush;
      }

      ctp_hdr_str(header);
      *ctp_out << header << "\n" << flush; 
   }
   
   // Create the CFP text file
   if(conf.output_flag(i_cfp).ival() >= 2) {

      // Create and open the CFP output file stream
      cfp_out = new ofstream;
      cfp_out->open(cfp_file);

      if(!cfp_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output CFP file \""
              << cfp_file << "\"\n\n";
         exit(1);
      }
      cfp_out->setf(ios::fixed);

      // List CFP file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output CFP file: "
              << cfp_file << "\n" << flush;
      }

      cfp_hdr_str(header);
      *cfp_out << header << "\n" << flush; 
   }
   
   // Create the COP text file
   if(conf.output_flag(i_cop).ival() >= 2) {

      // Create and open the COP output file stream
      cop_out = new ofstream;
      cop_out->open(cop_file);

      if(!cop_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output COP file \""
              << cop_file << "\"\n\n";
         exit(1);
      }
      cop_out->setf(ios::fixed);

      // List COP file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output COP file: "
              << cop_file << "\n" << flush;
      }

      cop_hdr_str(header);
      *cop_out << header << "\n" << flush; 
   }
   
   // Create the CTS text file
   if(conf.output_flag(i_cts).ival() >= 2) {

      // Create and open the CTS output file stream
      cts_out = new ofstream;
      cts_out->open(cts_file);

      if(!cts_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output CTS file \""
              << cts_file << "\"\n\n";
         exit(1);
      }
      cts_out->setf(ios::fixed);

      // List CTS file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output CTS file: "
              << cts_file << "\n" << flush;
      }
      
      cts_hdr_str(header);
      *cts_out << header << "\n" << flush; 
   }
   
   // Create the CNT text file
   if(conf.output_flag(i_cnt).ival() >= 2) {

      // Create and open the CNT output file stream
      cnt_out = new ofstream;
      cnt_out->open(cnt_file);

      if(!cnt_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output CNT file \""
              << cnt_file << "\"\n\n";
         exit(1);
      }
      cnt_out->setf(ios::fixed);

      // List CNT file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output CNT file: "
              << cnt_file << "\n" << flush;
      }
      
      cnt_hdr_str(header);
      *cnt_out << header << "\n" << flush; 
   }
   
   // Create the SL1L2 text file
   if(conf.output_flag(i_sl1l2).ival() >= 2) {

      // Create and open the SL1L2 output file stream
      sl1l2_out = new ofstream;
      sl1l2_out->open(sl1l2_file);

      if(!sl1l2_out) {
         cerr << "\n\nERROR: setup_txt_files() -> "
              << "unable to open output SL1L2 file \""
              << sl1l2_file << "\"\n\n";
         exit(1);
      }
      sl1l2_out->setf(ios::fixed);

      // List SL1L2 file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output SL1L2 file: "
              << sl1l2_file << "\n" << flush;
      }
      
      sl1l2_hdr_str(header);
      *sl1l2_out << header << "\n" << flush; 
   }   

   return;
}

////////////////////////////////////////////////////////////////////////

void setup_nc_file(unixtime valid_ut, int lead_sec) {
   int yr, mon, day, hr, min, sec;
   int l_hr, l_min, l_sec;
   int x, y;
   char attribute_str[max_str_len], time_str[max_str_len];
   char hostname_str[max_str_len];
   double lat, lon;
   float *lat_data, *lon_data;
   unixtime ut;

   // Create output NetCDF file name
   sec_to_hms(lead_sec, l_hr, l_min, l_sec);
   unix_to_mdyhms(valid_ut, mon, day, yr, hr, min, sec);
   sprintf(out_nc_file, "%s/grid_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_pairs.nc",
           out_dir, 
           l_hr, l_min, l_sec,
           yr, mon, day, hr, min, sec);

   // Create a new NetCDF file and open it
   nc_out = new NcFile(out_nc_file, NcFile::Replace);

   if(!nc_out || !nc_out->is_valid()) {
      cerr << "\n\nERROR: setup_nc_file() -> "
           << "trouble opening output NetCDF file "
           << out_nc_file << "\n\n" << flush;
      exit(1);
   }

   // List NetCDF file as it is being created
   if(verbosity > 1) {
      cout << "Creating Output NetCDF file: "
           << out_nc_file << "\n" << flush;
   }

   // Add global attributes
   nc_out->add_att("Conventions", "COARDS");

   ut = time(NULL);
   unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
   sprintf(time_str, "%.4i-%.2i-%.2i %.2i:%.2i:%.2i", yr, mon, day, hr, min, sec);

   gethostname(hostname_str, max_str_len);

   sprintf(attribute_str, "File %s generated %s UTC on host %s",
           out_nc_file, time_str, hostname_str);
   nc_out->add_att("FileOrigins", attribute_str);
   nc_out->add_att("Projection", grid.name());
   nc_out->add_att("Difference", "Forecast Value - Observation Value");

   // Define Dimensions
   lat_dim = nc_out->add_dim("lat", (long) grid.ny());
   lon_dim = nc_out->add_dim("lon", (long) grid.nx());

   // Define Variables
   lat_var = nc_out->add_var("lat", ncFloat, lat_dim, lon_dim);
   lon_var = nc_out->add_var("lon", ncFloat, lat_dim, lon_dim);

   // Add variable attributes
   lat_var->add_att("units", "degrees_north");
   lat_var->add_att("long_name", "Latitude");
   lat_var->add_att("_FillValue", bad_data_float);

   lon_var->add_att("units", "degrees_east");
   lon_var->add_att("long_name", "Longitude");
   lon_var->add_att("_FillValue", bad_data_float);

   // Allocate memory to store the Lat/Lon values for each grid point
   lat_data = lon_data = (float *) 0;
   lat_data = new float [grid.nx()*grid.ny()];
   lon_data = new float [grid.nx()*grid.ny()];

   for(x=0; x<grid.nx(); x++) {
      for(y=0; y<grid.ny(); y++) {

         // Lat/Lon values for each grid point
         grid.xy_to_latlon(x, y, lat, lon);

         lat_data[two_to_one(grid.nx(), grid.ny(), x, y)] = (float) lat;
         lon_data[two_to_one(grid.nx(), grid.ny(), x, y)] = (float) -1.0*lon;
      }
   }
   if( !lat_var->put(&lat_data[0], grid.ny(), grid.nx())
    || !lon_var->put(&lon_data[0], grid.ny(), grid.nx()) ) {

      cerr << "\n\nERROR: setup_nc_file() -> "
           << "error with the lat_var->put or lon_var->put\n\n" << flush;
      exit(1);
   }

   // Dealloate memory
   delete [] lat_data;
   lat_data = (float *) 0;
   delete [] lon_data;
   lon_data = (float *) 0;
   
   return;
}

////////////////////////////////////////////////////////////////////////

void write_txt_files(const WrfData &fcst_wd, const WrfData &obs_wd, int i_interp, int lev) {
   int i, j, mon, day, yr, hr, min, sec, status;
   int l_hr, l_min, l_sec;
   char line[max_line_len], fmt_str[max_line_len];
   char lead_str[max_str_len], valid_str[max_str_len], type_str[max_str_len];
   char thresh_str[max_str_len], var_str[max_num_vars][var_str_len];
   char interp_mthd_str[max_str_len], interp_pnts_str[max_str_len];
   unixtime ut;

   // Compute a string for the lead time and valid time
   ut = fcst_wd.get_valid_time();
   unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
   sec_to_hms(fcst_wd.get_lead_time(), l_hr, l_min, l_sec);
   sprintf(lead_str, "%.2i%.2i%.2i", l_hr, l_min, l_sec);
   sprintf(valid_str, "%.4i%.2i%.2i_%.2i%.2i%.2i", yr, mon, day, hr, min, sec);
   
   // Smoothing method and number of points
   strcpy(interp_mthd_str, mthd_str[interp_mthd[i_interp]]);
   sprintf(interp_pnts_str, "%i", interp_wdth[i_interp]*interp_wdth[i_interp]);

   // Dump out a VSDB record for each Contingency Table
   for(i=0; i<n_cts; i++) {

      // If the name of the CTS Info object is empty, continue
      if(!cts_info[i_interp][i].name ||
         strlen(cts_info[i_interp][i].name) <= 0) continue;

      // Type of thresholding performed
      if(cts_info[i_interp][i].cts_thresh_ind >= 0
      && cts_info[i_interp][i].cts_thresh_ind < n_thresh_type) {
         sprintf(thresh_str, "%s%.3f", 
            thresh_type_str[cts_info[i_interp][i].cts_thresh_ind], 
            cts_info[i_interp][i].cts_thresh);
      }
      else {
         cerr << "\n\nERROR: write_txt_files() -> "
              << "unexpected threshold indicator value of " 
              << cts_info[i_interp][i].cts_thresh_ind << "\".\n\n" << flush;
         exit(1);
      }

      // FHO
      // Dump out the FHO VSDB line: 
      //    TOTAL, F_RATE, H_RATE, O_RATE, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_fho).ival() != 0) {      
         
         sprintf(type_str, "FHO%s", thresh_str);

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
         dbl_to_str(cts_info[i_interp][i].cts.f_rate(status), var_str[1]);
         dbl_to_str(cts_info[i_interp][i].cts.h_rate(status), var_str[2]);
         dbl_to_str(cts_info[i_interp][i].cts.o_rate(status), var_str[3]);

         // Setup the format string
         fho_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,      // MET Version
            model_name,       // Model Name
            lead_str,         // Valid Time HHMMSS
            valid_str,        // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,        // VSDB code
            cts_info[i_interp][i].name, // Verification masking region
            type_str,         // Threshold applied
            field_name,       // Grib field name
            level_name,       // Grib level name
            var_str[0],       // Total number of points
            var_str[1],       // Forecast Rate = FY/N
            var_str[2],       // Hit Rate = FY_OY/N
            var_str[3],       // Observation = OY/N
            interp_mthd_str,  // Interpolation method 
            interp_pnts_str   // Interpolation points
         );
         *vsdb_out << line << "\n";
 
         if(conf.output_flag(i_fho).ival() >= 2) *fho_out << line << "\n";
      }

      // Contingency Table Counts
      // Dump out the CTC VSDB line:
      //    TOTAL, FY_OY, FY_ON, FN_OY, FN_ON, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_ctc).ival() != 0) {

         sprintf(type_str, "CTC%s", thresh_str);

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
         sprintf(var_str[1], "%i", cts_info[i_interp][i].cts.fy_oy());
         sprintf(var_str[2], "%i", cts_info[i_interp][i].cts.fy_on());
         sprintf(var_str[3], "%i", cts_info[i_interp][i].cts.fn_oy());
         sprintf(var_str[4], "%i", cts_info[i_interp][i].cts.fn_on());
      
         // Setup the format string
         ctc_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,      // MET Version
            model_name,       // Model Name
            lead_str,         // Valid Time HHMMSS
            valid_str,        // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,        // VSDB code
            cts_info[i_interp][i].name, // Verification masking region
            type_str,         // Threshold applied
            field_name,       // Grib field name
            level_name,       // Grib level name
            var_str[0],       // Total count
            var_str[1],       // fy_oy count
            var_str[2],       // fy_on count
            var_str[3],       // fn_oy count
            var_str[4],       // fn_on count
            interp_mthd_str,  // Interpolation method 
            interp_pnts_str   // Interpolation points
         );
         *vsdb_out << line << "\n";

         if(conf.output_flag(i_ctc).ival() >= 2) *ctc_out << line << "\n";
      }

      // Contingency Table Proportions
      // Dump out the CTP VSDB line:
      //    TOTAL, FY_OY_TP, FY_ON_TP, FN_OY_TP, FN_ON_TP, 
      //    FY_TP, FN_TP, OY_TP, ON_TP, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_ctp).ival() != 0) {
            
         sprintf(type_str, "CTP%s", thresh_str);

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
         dbl_to_str(cts_info[i_interp][i].cts.fy_oy_tp(status), var_str[1]);
         dbl_to_str(cts_info[i_interp][i].cts.fy_on_tp(status), var_str[2]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_oy_tp(status), var_str[3]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_on_tp(status), var_str[4]);
         dbl_to_str(cts_info[i_interp][i].cts.fy_tp(status), var_str[5]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_tp(status), var_str[6]);
         dbl_to_str(cts_info[i_interp][i].cts.oy_tp(status), var_str[7]);
         dbl_to_str(cts_info[i_interp][i].cts.on_tp(status), var_str[8]);

         // Setup the format string
         ctp_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,      // MET Version
            model_name,       // Model Name
            lead_str,         // Valid Time HHMMSS
            valid_str,        // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,        // VSDB code
            cts_info[i_interp][i].name, // Verification masking region
            type_str,         // Threshold applied
            field_name,       // Grib field name
            level_name,       // Grib level name
            var_str[0],       // Total count
            var_str[1],       // fy_oy_tp
            var_str[2],       // fy_on_tp
            var_str[3],       // fn_oy_tp
            var_str[4],       // fn_on_tp
            var_str[5],       // fy_tp
            var_str[6],       // fn_tp
            var_str[7],       // oy_tp
            var_str[8],       // on_tp
            interp_mthd_str,  // Interpolation method 
            interp_pnts_str   // Interpolation points
         );
         *vsdb_out << line << "\n";

         if(conf.output_flag(i_ctp).ival() >= 2) *ctp_out << line << "\n";
      }

      // Contingency Table Forecast Proportions
      // Dump out the CFP VSDB line:
      //    TOTAL, FY_OY_FP, FY_ON_FP, FN_OY_FP, FN_ON_FP, FY, FN, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_cfp).ival() != 0) {
      
         sprintf(type_str, "CFP%s", thresh_str);

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
         dbl_to_str(cts_info[i_interp][i].cts.fy_oy_fp(status), var_str[1]);
         dbl_to_str(cts_info[i_interp][i].cts.fy_on_fp(status), var_str[2]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_oy_fp(status), var_str[3]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_on_fp(status), var_str[4]);
         sprintf(var_str[5], "%i", cts_info[i_interp][i].cts.fy());
         sprintf(var_str[6], "%i", cts_info[i_interp][i].cts.fn());

         // Setup the format string
         cfp_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,      // MET Version
            model_name,       // Model Name
            lead_str,         // Valid Time HHMMSS
            valid_str,        // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,        // VSDB code
            cts_info[i_interp][i].name, // Verification masking region
            type_str,         // Threshold applied
            field_name,       // Grib field name
            level_name,       // Grib level name
            var_str[0],       // Total count
            var_str[1],       // fy_oy_fp
            var_str[2],       // fy_on_fp
            var_str[3],       // fn_oy_fp
            var_str[4],       // fn_on_fp
            var_str[5],       // fy count
            var_str[6],       // fn count
            interp_mthd_str,  // Interpolation method 
            interp_pnts_str   // Interpolation points
         );
         *vsdb_out << line << "\n";

         if(conf.output_flag(i_cfp).ival() >= 2) *cfp_out << line << "\n";
      }
      
      // Contingency Table Observation Proportions
      // Dump out the COP VSDB line:
      //    TOTAL, FY_OY_OP, FY_ON_OP, FN_OY_OP, FN_ON_OP, OY, ON, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_cop).ival() != 0) {
      
         sprintf(type_str, "COP%s", thresh_str);

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
         dbl_to_str(cts_info[i_interp][i].cts.fy_oy_op(status), var_str[1]);
         dbl_to_str(cts_info[i_interp][i].cts.fy_on_op(status), var_str[2]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_oy_op(status), var_str[3]);
         dbl_to_str(cts_info[i_interp][i].cts.fn_on_op(status), var_str[4]);
         sprintf(var_str[5], "%i", cts_info[i_interp][i].cts.oy());
         sprintf(var_str[6], "%i", cts_info[i_interp][i].cts.on());

         // Setup the format string
         cop_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,          // MET Version
            model_name,           // Model Name
            lead_str,             // Valid Time HHMMSS
            valid_str,            // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,            // VSDB code
            cts_info[i_interp][i].name,     // Verification masking region
            type_str,             // Threshold applied
            field_name,           // Grib field name
            level_name,           // Grib level name
            var_str[0],           // Total count
            var_str[1],           // fy_oy_op
            var_str[2],           // fy_on_op
            var_str[3],           // fn_oy_op
            var_str[4],           // fn_on_op
            var_str[5],           // oy count
            var_str[6],           // on count
            interp_mthd_str,      // Interpolation method 
            interp_pnts_str       // Interpolation points
         );
         *vsdb_out << line << "\n";

         if(conf.output_flag(i_cop).ival() >= 2) *cop_out << line << "\n";
      }
      
      // Contingency Table Stats
      // Dump out the CTS VSDB line: 
      //    TOTAL, 
      //    BASER, BASER_CL, BASER_CU,
      //    FMEAN, FMEAN_CL, FMEAN_CU,
      //    ACC, ACC_CL, ACC_CU,
      //    FBIAS,
      //    PODY, PODY_CL, PODY_CU,
      //    PODN, PODN_CL, PODN_CU,
      //    POFD, POFD_CL, POFD_CU,
      //    FAR, FAR_CL, FAR_CU,
      //    CSI, CSI_CL, CSI_CU,
      //    GSS,
      //    HK, HK_CL, HK_CU,
      //    HSS,
      //    ODDS, ODDS_CL, ODDS_CU
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_cts).ival() != 0) {
      
         for(j=0; j<cts_info[i_interp][i].n_alpha; j++) {

            sprintf(type_str, "CTS%s/%.3f", thresh_str, cts_info[i_interp][i].alpha[j]);

            // Convert int and double values to strings
            sprintf(var_str[0], "%i", cts_info[i_interp][i].cts.n());
            dbl_to_str(cts_info[i_interp][i].baser.v, var_str[1]);
            dbl_to_str(cts_info[i_interp][i].baser.v_cl[j], var_str[2]);
            dbl_to_str(cts_info[i_interp][i].baser.v_cu[j], var_str[3]);
            dbl_to_str(cts_info[i_interp][i].fmean.v, var_str[4]);
            dbl_to_str(cts_info[i_interp][i].fmean.v_cl[j], var_str[5]);
            dbl_to_str(cts_info[i_interp][i].fmean.v_cu[j], var_str[6]);
            dbl_to_str(cts_info[i_interp][i].acc.v, var_str[7]);
            dbl_to_str(cts_info[i_interp][i].acc.v_cl[j], var_str[8]);
            dbl_to_str(cts_info[i_interp][i].acc.v_cu[j], var_str[9]);
            dbl_to_str(cts_info[i_interp][i].fbias.v, var_str[10]);
            dbl_to_str(cts_info[i_interp][i].pody.v, var_str[11]);
            dbl_to_str(cts_info[i_interp][i].pody.v_cl[j], var_str[12]);
            dbl_to_str(cts_info[i_interp][i].pody.v_cu[j], var_str[13]);
            dbl_to_str(cts_info[i_interp][i].podn.v, var_str[14]);
            dbl_to_str(cts_info[i_interp][i].podn.v_cl[j], var_str[15]);
            dbl_to_str(cts_info[i_interp][i].podn.v_cu[j], var_str[16]);
            dbl_to_str(cts_info[i_interp][i].pofd.v, var_str[17]);
            dbl_to_str(cts_info[i_interp][i].pofd.v_cl[j], var_str[18]);
            dbl_to_str(cts_info[i_interp][i].pofd.v_cu[j], var_str[19]);
            dbl_to_str(cts_info[i_interp][i].far.v, var_str[20]);
            dbl_to_str(cts_info[i_interp][i].far.v_cl[j], var_str[21]);
            dbl_to_str(cts_info[i_interp][i].far.v_cu[j], var_str[22]);
            dbl_to_str(cts_info[i_interp][i].csi.v, var_str[23]);
            dbl_to_str(cts_info[i_interp][i].csi.v_cl[j], var_str[24]);
            dbl_to_str(cts_info[i_interp][i].csi.v_cu[j], var_str[25]);
            dbl_to_str(cts_info[i_interp][i].gss.v, var_str[26]);
            dbl_to_str(cts_info[i_interp][i].hk.v, var_str[27]);
            dbl_to_str(cts_info[i_interp][i].hk.v_cl[j], var_str[28]);
            dbl_to_str(cts_info[i_interp][i].hk.v_cu[j], var_str[29]);
            dbl_to_str(cts_info[i_interp][i].hss.v, var_str[30]);
            dbl_to_str(cts_info[i_interp][i].odds.v, var_str[31]);
            dbl_to_str(cts_info[i_interp][i].odds.v_cl[j], var_str[32]);
            dbl_to_str(cts_info[i_interp][i].odds.v_cu[j], var_str[33]);

            // Setup the format string
            cts_fmt_str(fmt_str);

            // Write the line
            sprintf(line, fmt_str,
               met_version,                           // MET Version
               model_name,                            // Model Name
               lead_str,                              // Valid Time HHMMSS
               valid_str,                             // Valid Time YYYYMMDD_HHMMSS
               vsdb_code,                             // VSDB code
               cts_info[i_interp][i].name,            // Verification masking region
               type_str,                              // Threshold applied
               field_name,                            // Grib field name
               level_name,                            // Grib level name
               var_str[0],                            // Total count
               var_str[1],   var_str[2],  var_str[3], // Base Rate (oy_tp) CI
               var_str[4],   var_str[5],  var_str[6], // Forecast Mean (fy_tp) CI
               var_str[7],   var_str[8],  var_str[9], // Accuracy CI
               var_str[10],                           // Frequency Bias
               var_str[11], var_str[12], var_str[13], // Probability of detecting yes CI
               var_str[14], var_str[15], var_str[16], // Probability of detecting no CI
               var_str[17], var_str[18], var_str[19], // Probability of false detection CI
               var_str[20], var_str[21], var_str[22], // False alarm rate CI
               var_str[23], var_str[24], var_str[25], // Critical Success Index (Threat Score) CI
               var_str[26],                           // Gilbert Skill Score (Equitable Threat Score)
               var_str[27], var_str[28], var_str[29], // Hanssen and Kuipers Discriminant (TSS) CI
               var_str[30],                           // Heidke Skill Score
               var_str[31], var_str[32], var_str[33], // Odds Ratio CI
               interp_mthd_str,                       // Interpolation method 
               interp_pnts_str                        // Interpolation points
            );
            *vsdb_out << line << "\n";

            if(conf.output_flag(i_cts).ival() >= 2) *cts_out << line << "\n";
         } // end for j
      } // end if
   } // end for i

   // Dump out a VSDB record for each Continuous Statistics object
   for(i=0; i<n_cnt; i++) {
      
      // Continuous Variable Stats
      // Dump out the CNT VSDB line: 
      //    TOTAL, 
      //    FBAR, FBAR_CL, FBAR_CU,
      //    FSTDEV, FSTDEV_CL, FSTDEV_CU,
      //    OBAR, OBAR_CL, OBAR_CU,
      //    OSTDEV, OSTDEV_CL, OSTDEV_CU,
      //    PR_CORR, PR_CORR_CL, PR_CORR_CU,
      //    SP_CORR, KT_CORR,
      //    RANKS, FRANK_TIES, ORANK_TIES,
      //    ME, ME_CL, ME_CU,
      //    ESTDEV, ESTDEV_CL, ESTDEV_CU,
      //    BIAS, MBIAS, MAE, MSE, BCMSE, RMSE,
      //    E10, E25, E50, E75, E90
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_cnt).ival() != 0) {
      
         for(j=0; j<cnt_info[i_interp][i].n_alpha; j++) {
         
            sprintf(type_str, "CNT/%.3f", cnt_info[i_interp][i].alpha[j]);

            // Convert int and double values to strings
            sprintf(var_str[0], "%i", cnt_info[i_interp][i].n);
            dbl_to_str(cnt_info[i_interp][i].fbar.v, var_str[1]);
            dbl_to_str(cnt_info[i_interp][i].fbar.v_cl[j], var_str[2]);
            dbl_to_str(cnt_info[i_interp][i].fbar.v_cu[j], var_str[3]);
            dbl_to_str(cnt_info[i_interp][i].fstdev.v, var_str[4]);
            dbl_to_str(cnt_info[i_interp][i].fstdev.v_cl[j], var_str[5]);
            dbl_to_str(cnt_info[i_interp][i].fstdev.v_cu[j], var_str[6]);
            dbl_to_str(cnt_info[i_interp][i].obar.v, var_str[7]);
            dbl_to_str(cnt_info[i_interp][i].obar.v_cl[j], var_str[8]);
            dbl_to_str(cnt_info[i_interp][i].obar.v_cu[j], var_str[9]);
            dbl_to_str(cnt_info[i_interp][i].ostdev.v, var_str[10]);
            dbl_to_str(cnt_info[i_interp][i].ostdev.v_cl[j], var_str[11]);
            dbl_to_str(cnt_info[i_interp][i].ostdev.v_cu[j], var_str[12]);
            dbl_to_str(cnt_info[i_interp][i].pr_corr.v, var_str[13]);
            dbl_to_str(cnt_info[i_interp][i].pr_corr.v_cl[j], var_str[14]);
            dbl_to_str(cnt_info[i_interp][i].pr_corr.v_cu[j], var_str[15]);
            dbl_to_str(cnt_info[i_interp][i].sp_corr, var_str[16]);
            dbl_to_str(cnt_info[i_interp][i].kt_corr, var_str[17]);
            sprintf(var_str[18], "%i", cnt_info[i_interp][i].n_ranks);
            sprintf(var_str[19], "%i", cnt_info[i_interp][i].frank_ties);
            sprintf(var_str[20], "%i", cnt_info[i_interp][i].orank_ties);
            dbl_to_str(cnt_info[i_interp][i].me.v, var_str[21]);
            dbl_to_str(cnt_info[i_interp][i].me.v_cl[j], var_str[22]);
            dbl_to_str(cnt_info[i_interp][i].me.v_cu[j], var_str[23]);
            dbl_to_str(cnt_info[i_interp][i].estdev.v, var_str[24]);
            dbl_to_str(cnt_info[i_interp][i].estdev.v_cl[j], var_str[25]);
            dbl_to_str(cnt_info[i_interp][i].estdev.v_cu[j], var_str[26]);
            dbl_to_str(cnt_info[i_interp][i].bias, var_str[27]);
            dbl_to_str(cnt_info[i_interp][i].mbias, var_str[28]);
            dbl_to_str(cnt_info[i_interp][i].mae, var_str[29]);
            dbl_to_str(cnt_info[i_interp][i].mse, var_str[30]);
            dbl_to_str(cnt_info[i_interp][i].bcmse, var_str[31]);
            dbl_to_str(cnt_info[i_interp][i].rmse, var_str[32]);
            dbl_to_str(cnt_info[i_interp][i].e10, var_str[33]);
            dbl_to_str(cnt_info[i_interp][i].e25, var_str[34]);
            dbl_to_str(cnt_info[i_interp][i].e50, var_str[35]);
            dbl_to_str(cnt_info[i_interp][i].e75, var_str[36]);
            dbl_to_str(cnt_info[i_interp][i].e90, var_str[37]);

            // Setup the format string
            cnt_fmt_str(fmt_str);

            // Write the line
            sprintf(line, fmt_str,
               met_version,                           // MET Version
               model_name,                            // Model Name
               lead_str,                              // Valid Time HHMMSS
               valid_str,                             // Valid Time YYYYMMDD_HHMMSS
               vsdb_code,                             // VSDB code
               cts_info[i_interp][i].name,            // Verification masking region
               type_str,                              // Threshold applied
               field_name,                            // Grib field name
               level_name,                            // Grib level name
               var_str[0],                            // Total Number of Grid Points
               var_str[1],   var_str[2],  var_str[3], // Forecast Mean CI
               var_str[4],   var_str[5],  var_str[6], // Forecast Standard Deviation CI
               var_str[7],   var_str[8],  var_str[9], // Observation Mean CI
               var_str[10], var_str[11], var_str[12], // Observation Standard Deviation CI
               var_str[13], var_str[14], var_str[15], // Pearson's Correlation Coefficient CI
               var_str[16],                           // Spearman's Rank Correlation Coefficient
               var_str[17],                           // Kendall Tau Rank Correlation Coefficient
               var_str[18],                           // Number of ranks compared
               var_str[19],                           // Number of tied ranks in the forecast field    
               var_str[20],                           // Number of tied ranks in the observation field
               var_str[21], var_str[22], var_str[23], // Mean Error CI
               var_str[24], var_str[25], var_str[26], // Error Standard Deviation CI
               var_str[27],                           // Bias
               var_str[28],                           // Multiplicative Bias
               var_str[29],                           // Mean Absolute Error         
               var_str[30],                           // Mean Squared Error
               var_str[31],                           // Bias-Corrected Mean Squared Error
               var_str[32],                           // Root Mean Squared Error
               var_str[33],                           // 10th Percentile of the Error
               var_str[34],                           // 25th Percentile of the Error
               var_str[35],                           // 50th Percentile of the Error
               var_str[36],                           // 75th Percentile of the Error
               var_str[37],                           // 90th Percentile of the Error
               interp_mthd_str,                       // Interpolation method 
               interp_pnts_str                        // Interpolation points
            );
            *vsdb_out << line << "\n";

            if(conf.output_flag(i_cnt).ival() >= 2) *cnt_out << line << "\n";
         } // end for j
      } // end if
      
      // Scalar L1L2 Line Type (SL1L2)
      // Dump out the SL1L2 VSDB line:
      //    TOTAL, FBAR, OBAR, FOBAR, FFBAR, OOBAR, 
      //    INTERP_MTHD, INTERP_PNTS
      if(conf.output_flag(i_sl1l2).ival() != 0) {
      
         strcpy(type_str, "SL1L2");

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cnt_info[i_interp][i].n);
         dbl_to_str(cnt_info[i_interp][i].fbar.v, var_str[1]);
         dbl_to_str(cnt_info[i_interp][i].obar.v, var_str[2]);
         dbl_to_str(cnt_info[i_interp][i].fobar, var_str[3]);
         dbl_to_str(cnt_info[i_interp][i].ffbar, var_str[4]);
         dbl_to_str(cnt_info[i_interp][i].oobar, var_str[5]);

         // Setup the format string
         sl1l2_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,      // MET Version
            model_name,       // Model Name
            lead_str,         // Valid Time HHMMSS
            valid_str,        // Valid Time YYYYMMDD_HHMMSS
            vsdb_code,        // VSDB code
            cts_info[i_interp][i].name, // Verification masking region
            type_str,         // Threshold applied
            field_name,       // Grib field name
            level_name,       // Grib level name
            var_str[0],       // Total count
            var_str[1],       // fbar
            var_str[2],       // obar
            var_str[3],       // fobar
            var_str[4],       // ffbar
            var_str[5],       // oobar
            interp_mthd_str,  // Interpolation method 
            interp_pnts_str   // Interpolation points
         );
         *vsdb_out << line << "\n";

         if(conf.output_flag(i_sl1l2).ival() >= 2) *sl1l2_out << line << "\n";
      }
   } // end for i
   
   return;
}

////////////////////////////////////////////////////////////////////////

void close_txt_files() {

   // Close the output VSDB file
   if(vsdb_out) {

      // List VSDB file after it is finished
      if(verbosity > 0) {
         cout << "Output VSDB file: "
              << vsdb_file << "\n" << flush;
      }
   
      // Close the output VSDB file
      vsdb_out->close(); 
      delete vsdb_out; 
      vsdb_out = (ofstream *) 0;
   }
   
   // Close the output FHO txt file
   if(conf.output_flag(i_fho).ival() >= 2 && fho_out) {

      // List FHO file after it is finished
      if(verbosity > 0) {
         cout << "Output FHO file: "
              << fho_file << "\n" << flush;
      }
   
      // Close the output FHO file
      fho_out->close(); 
      delete fho_out; 
      fho_out = (ofstream *) 0;
   }
   
   // Close the output CTC txt file
   if(conf.output_flag(i_ctc).ival() >= 2 && ctc_out) {

      // List CTC file after it is finished
      if(verbosity > 0) {
         cout << "Output CTC file: "
              << ctc_file << "\n" << flush;
      }
   
      // Close the output CTC file
      ctc_out->close(); 
      delete ctc_out; 
      ctc_out = (ofstream *) 0;
   }

   // Close the output CTP txt file
   if(conf.output_flag(i_ctp).ival() >= 2 && ctp_out) {

      // List CTP file after it is finished
      if(verbosity > 0) {
         cout << "Output CTP file: "
              << ctp_file << "\n" << flush;
      }
   
      // Close the output CTP file
      ctp_out->close(); 
      delete ctp_out; 
      ctp_out = (ofstream *) 0;
   }
   
   // Close the output CFP txt file
   if(conf.output_flag(i_cfp).ival() >= 2 && cfp_out) {

      // List CFP file after it is finished
      if(verbosity > 0) {
         cout << "Output CFP file: "
              << cfp_file << "\n" << flush;
      }
   
      // Close the output CFP file
      cfp_out->close(); 
      delete cfp_out; 
      cfp_out = (ofstream *) 0;
   }
   
   // Close the output COP txt file
   if(conf.output_flag(i_cop).ival() >= 2 && cop_out) {

      // List COP file after it is finished
      if(verbosity > 0) {
         cout << "Output COP file: "
              << cop_file << "\n" << flush;
      }
   
      // Close the output COP file
      cop_out->close(); 
      delete cop_out; 
      cop_out = (ofstream *) 0;
   }
      
   // Close the output CTS txt file
   if(conf.output_flag(i_cts).ival() >= 2 && cts_out) {

      // List CTS file after it is finished
      if(verbosity > 0) {
         cout << "Output CTS file: "
              << cts_file << "\n" << flush;
      }
   
      // Close the output CTS file
      cts_out->close(); 
      delete cts_out; 
      cts_out = (ofstream *) 0;
   }
   
   // Close the output CNT txt file
   if(conf.output_flag(i_cnt).ival() >= 2 && cnt_out) {

      // List CNT file after it is finished
      if(verbosity > 0) {
         cout << "Output CNT file: "
              << cnt_file << "\n" << flush;
      }
   
      // Close the output CNT file
      cnt_out->close(); 
      delete cnt_out; 
      cnt_out = (ofstream *) 0;
   }

   // Close the output SL1L2 txt file
   if(conf.output_flag(i_sl1l2).ival() >= 2 && sl1l2_out) {

      // List SL1L2 file after it is finished
      if(verbosity > 0) {
         cout << "Output SL1L2 file: "
              << sl1l2_file << "\n" << flush;
      }
   
      // Close the output SL1L2 file
      sl1l2_out->close(); 
      delete sl1l2_out; 
      sl1l2_out = (ofstream *) 0;
   }
      
   return;
}

////////////////////////////////////////////////////////////////////////

void compute_cts(const WrfData &fcst_wd, const WrfData &obs_wd, int i_interp, int i_gc) {
   int  i, j, k, n, x, y, status;
   WrfData fwd, owd;
   
   if(verbosity > 1)
      cout << "Computing Scores for [Field, Threshold, "
           << "Vx Region] =\n" << flush;

   // Initialize the CTS Info object name to an empty string
   for(i=0; i<n_cts; i++) cts_info[i_interp][i].set_name("");

   // Score the fields for each combination of threshold and 
   // verification region specified
   for(i=0; i<gc_ti[i_gc].n_thresh; i++) {
      
      // Initialize the two fields
      fwd = fcst_wd;
      owd = obs_wd;
      
      // Threshold the raw forecast and observation fields using
      // the type of thresholding indicated
      fwd.threshold_double(gc_ti[i_gc].thresh[i], gc_ti[i_gc].thresh_ind[i]);
      owd.threshold_double(gc_ti[i_gc].thresh[i], gc_ti[i_gc].thresh_ind[i]);

      for(j=0; j<n_masks; j++) {
       
         n = i*n_masks + j;
            
         if(verbosity > 1) 
            cout << "[" << field_name
                 << ", " << gc_ti[i_gc].thresh[i]
                 << ", " << mask_names[j] << "]\n" << flush; 
      
         cts_info[i_interp][n].cts_thresh = gc_ti[i_gc].thresh[i];
         cts_info[i_interp][n].cts_thresh_ind = gc_ti[i_gc].thresh_ind[i];
         cts_info[i_interp][n].set_name(mask_names[j]);
         cts_info[i_interp][n].cts.zero_out();

         for(x=0; x<fwd.get_nx(); x++) {
            for(y=0; y<fwd.get_ny(); y++) {

               // Skip any grid points containing bad data values for
               // either of the raw fields or where the verification
               // masking region is turned off
               if(fcst_wd.is_bad_data(x, y) || 
                  obs_wd.is_bad_data(x, y) ||
                  !mask_wd[j].s_is_on(x, y) ) continue;

               else if( fwd.s_is_on(x, y) &&  owd.s_is_on(x, y)) cts_info[i_interp][n].cts.inc_fy_oy();
               else if( fwd.s_is_on(x, y) && !owd.s_is_on(x, y)) cts_info[i_interp][n].cts.inc_fy_on();
               else if(!fwd.s_is_on(x, y) &&  owd.s_is_on(x, y)) cts_info[i_interp][n].cts.inc_fn_oy();
               else if(!fwd.s_is_on(x, y) && !owd.s_is_on(x, y)) cts_info[i_interp][n].cts.inc_fn_on();
            }
         }

         // Compute contingency table statistics
         cts_info[i_interp][n].baser.v = cts_info[i_interp][n].cts.oy_tp(status);
         cts_info[i_interp][n].fmean.v = cts_info[i_interp][n].cts.fy_tp(status);
         cts_info[i_interp][n].acc.v   = cts_info[i_interp][n].cts.accuracy(status);
         cts_info[i_interp][n].fbias.v = cts_info[i_interp][n].cts.fbias(status);
         cts_info[i_interp][n].pody.v  = cts_info[i_interp][n].cts.pod_yes(status);
         cts_info[i_interp][n].podn.v  = cts_info[i_interp][n].cts.pod_no(status);
         cts_info[i_interp][n].pofd.v  = cts_info[i_interp][n].cts.pofd(status);
         cts_info[i_interp][n].far.v   = cts_info[i_interp][n].cts.far(status);
         cts_info[i_interp][n].csi.v   = cts_info[i_interp][n].cts.csi(status);
         cts_info[i_interp][n].gss.v   = cts_info[i_interp][n].cts.gss(status);
         cts_info[i_interp][n].hk.v    = cts_info[i_interp][n].cts.hk(status);
         cts_info[i_interp][n].hss.v   = cts_info[i_interp][n].cts.hss(status);
         cts_info[i_interp][n].odds.v  = cts_info[i_interp][n].cts.odds(status);
      
         // For each choice of alpha, compute the corresponding 
         // confidence interval
         for(k=0; k<cts_info[i_interp][n].n_alpha; k++) {

            compute_proportion_ci(cts_info[i_interp][n].baser.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].baser.v_cl[k], cts_info[i_interp][n].baser.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].fmean.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].fmean.v_cl[k], cts_info[i_interp][n].fmean.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].acc.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].acc.v_cl[k], cts_info[i_interp][n].acc.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].pody.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].pody.v_cl[k], cts_info[i_interp][n].pody.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].podn.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].podn.v_cl[k], cts_info[i_interp][n].podn.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].pofd.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].pofd.v_cl[k], cts_info[i_interp][n].pofd.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].far.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].far.v_cl[k], cts_info[i_interp][n].far.v_cu[k]);
            compute_proportion_ci(cts_info[i_interp][n].csi.v, cts_info[i_interp][n].cts.n(), cts_info[i_interp][n].alpha[k], 
                                  cts_info[i_interp][n].csi.v_cl[k], cts_info[i_interp][n].csi.v_cu[k]);

            // Compute a confidence interval for Hanssen and Kuipers discriminant
            compute_hk_ci(cts_info[i_interp][n].hk.v, cts_info[i_interp][n].alpha[k], 
                          cts_info[i_interp][n].cts.fy_oy(), cts_info[i_interp][n].cts.fy_on(), 
                          cts_info[i_interp][n].cts.fn_oy(), cts_info[i_interp][n].cts.fn_on(), 
                          cts_info[i_interp][n].hk.v_cl[k], cts_info[i_interp][n].hk.v_cu[k]);

            // Compute a confidence interval for the odds ratio
            compute_woolf_ci(cts_info[i_interp][n].odds.v, cts_info[i_interp][n].alpha[k], 
                             cts_info[i_interp][n].cts.fy_oy(), cts_info[i_interp][n].cts.fy_on(), 
                             cts_info[i_interp][n].cts.fn_oy(), cts_info[i_interp][n].cts.fn_on(), 
                             cts_info[i_interp][n].odds.v_cl[k], cts_info[i_interp][n].odds.v_cu[k]);
         } // end for k
      } // end for j
   } // end for i

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Compute continuous statistics for the fields provided.
// "gc" indicates the grib code being verified.
//
////////////////////////////////////////////////////////////////////////

void compute_cnt(const WrfData &fcst_wd, const WrfData &obs_wd, int i_interp, int gc) {
   int i, j, k, x, y, count;
   int n_zero_zero, n_fcst_rank, n_obs_rank, n_fcst_rank_ties, n_obs_rank_ties;
   int concordant, discordant, extra_f, extra_o;
   double f, o, f_sum, o_sum, ff_sum, oo_sum, fo_sum;
   double err, err_sum, abs_err_sum, err_sq_sum, den;
   double cv_normal_l, cv_normal_u, cv_chi2_l, cv_chi2_u, v, cl, cu;
   double *err_array, *fcst_rank, *obs_rank;
   WrfData fcst_wd_mask_bad, obs_wd_mask_bad;
   WrfData fcst_wd_mask, obs_wd_mask;
   WrfData fcst_wd_rank, obs_wd_rank;
   
   if(verbosity > 1)
      cout << "Computing Continuous Scores for [Field, Vx Region, Alpha] =\n" 
           << flush;

   // Use the bad data in the forecast field to mask the obs field
   // and vice versa, for use in computing the ranks
   fcst_wd_mask_bad = fcst_wd;
   obs_wd_mask_bad = obs_wd;
   mask_bad_data(fcst_wd_mask_bad, obs_wd);
   mask_bad_data(obs_wd_mask_bad, fcst_wd);

   // Allocate space to store the differences for use in computing the 
   // percentiles and for storing the data ranks
   err_array = new double [fcst_wd.get_nx()*fcst_wd.get_ny()];
   fcst_rank = new double [fcst_wd.get_nx()*fcst_wd.get_ny()];
   obs_rank  = new double [fcst_wd.get_nx()*fcst_wd.get_ny()];
   
   // Compute the CNT scores for each masking region and choice
   // of alpha
   for(i=0; i<n_cnt; i++) {
      
      // Store the name of the mask in the CNTInfo object
      cnt_info[i_interp][i].set_name(mask_names[i]);

      // Compute the difference field and sums of the differences
      count = 0;
      f_sum = o_sum = ff_sum = oo_sum = fo_sum = 0.0;
      err_sum = abs_err_sum = err_sq_sum = 0.0;
      for(x=0; x<fcst_wd.get_nx(); x++) {
         for(y=0; y<fcst_wd.get_ny(); y++) {

            // Skip any grid points containing bad data values for
            // either of the raw fields or where the verification
            // masking region is turned off
            if(fcst_wd.is_bad_data(x, y) || obs_wd.is_bad_data(x, y) ||
               !mask_wd[i].s_is_on(x, y) ) continue;
            
            f = fcst_wd.get_xy_double(x, y);
            o = obs_wd.get_xy_double(x, y);
 
            // Compute the error
            err = f-o;
            err_array[count] = err;

            f_sum       += f;
            o_sum       += o;
            ff_sum      += f*f;
            oo_sum      += o*o;
            fo_sum      += f*o;
            err_sum     += err;
            abs_err_sum += abs(err);
            err_sq_sum  += err*err;
            count       += 1;
         }
      }
      
      if(count == 0) {
         cerr << "\n\nERROR: compute_cnt() -> "
              << "count is zero!\n\n" << flush;
         exit(1);
      }

      // Store the sample size
      cnt_info[i_interp][i].n = count;
      
      // Compute the forecast mean
      cnt_info[i_interp][i].fbar.v = f_sum/count;

      // Compute the forecast standard deviation
      cnt_info[i_interp][i].fstdev.v = compute_stdev(f_sum, ff_sum, count);

      // Compute the observation mean
      cnt_info[i_interp][i].obar.v = o_sum/count;

      // Compute the observation standard deviation
      cnt_info[i_interp][i].ostdev.v = compute_stdev(o_sum, oo_sum, count);

      // Compute the f*o mean
      cnt_info[i_interp][i].fobar = fo_sum/count;
      
      // Compute the forecast squared mean
      cnt_info[i_interp][i].ffbar = ff_sum/count;
      
      // Compute the observation squared mean
      cnt_info[i_interp][i].oobar = oo_sum/count;
      
      // Compute the bias
      cnt_info[i_interp][i].bias = cnt_info[i_interp][i].fbar.v - cnt_info[i_interp][i].obar.v;
      
      // Compute multiplicative bias
      if(cnt_info[i_interp][i].obar.v == 0) 
         cnt_info[i_interp][i].mbias = bad_data_double;
      else
         cnt_info[i_interp][i].mbias = cnt_info[i_interp][i].fbar.v/cnt_info[i_interp][i].obar.v;

      // Compute the correlation coefficient
      den = sqrt((count*ff_sum - f_sum*f_sum)*(count*oo_sum - o_sum*o_sum));
      if(den == 0.0) cnt_info[i_interp][i].pr_corr.v = bad_data_double; 
      else           cnt_info[i_interp][i].pr_corr.v = ((count*fo_sum) - (f_sum*o_sum))/den;

      //
      // Set up field for computing the Kendall Tau and Spearman Rank 
      // correlation coefficients
      //
      
      // Initialize the fcst and obs mask field to the mask bad data fields
      fcst_wd_mask = fcst_wd_mask_bad;
      obs_wd_mask = obs_wd_mask_bad;
      
      // Apply the current mask field to raw fcst and obs fields
      apply_mask(fcst_wd_mask, mask_wd[i]);
      apply_mask(obs_wd_mask, mask_wd[i]);
      
      // If verifying precipitation, mask out the (0, 0) cases
      // Apply (0,0) mask for:
      // - precipitation rate
      // - thunderstorm probability
      // - total precipitation
      // - large scale precipitation
      // - convective precipitation
      if(gc == prate_grib_code || gc == tstm_grib_code  || 
         gc == apcp_grib_code  || gc == ncpcp_grib_code || 
         gc == acpcp_grib_code) {
         n_zero_zero = mask_double_double(fcst_wd_mask, obs_wd_mask, 0.0);
      }
      else {
         n_zero_zero = 0;
      }
      
      // Compute the rank of the remaining raw data values 
      // in the fcst and obs fields
      n_fcst_rank = compute_rank(fcst_wd_mask, fcst_wd_rank, 
                                 fcst_rank, n_fcst_rank_ties);
      n_obs_rank  = compute_rank(obs_wd_mask, obs_wd_rank, 
                                 obs_rank, n_obs_rank_ties);
      
      if(n_fcst_rank != n_obs_rank) {
         cerr << "\n\nERROR: compute_cnt() -> "
              << "the number of ranked data points in the forecast field ("
              << n_fcst_rank << ") does not equal the number of ranked data points "
              << "in the obs field (" << n_obs_rank << ")\n\n" << flush;
         exit(1);      
      }
      cnt_info[i_interp][i].n_ranks = n_fcst_rank;
      cnt_info[i_interp][i].frank_ties = n_fcst_rank_ties;
      cnt_info[i_interp][i].orank_ties = n_obs_rank_ties;
      
      // Compute the sums for the ranks for use in computing Spearman's
      // Rank correlation coefficient
      count = 0;
      f_sum = o_sum = ff_sum = oo_sum = fo_sum = 0.0;
      for(x=0; x<fcst_wd_rank.get_nx(); x++) {
         for(y=0; y<fcst_wd_rank.get_ny(); y++) {

            // Skip any grid points containing bad data values for
            // either of the raw fields or where the verification
            // masking region is turned off
            if(fcst_wd_rank.is_bad_data(x, y) || 
               obs_wd_rank.is_bad_data(x, y)) continue;
            
            f = fcst_wd_rank.get_xy_double(x, y);
            o = obs_wd_rank.get_xy_double(x, y);

            f_sum  += f;
            o_sum  += o;
            ff_sum += f*f;
            oo_sum += o*o;
            fo_sum += f*o;
            count  += 1;
         }
      }
      
      // Compute Spearman's Rank correlation coefficient
      den = sqrt((count*ff_sum - f_sum*f_sum)*(count*oo_sum - o_sum*o_sum));
      if(den == 0.0) cnt_info[i_interp][i].sp_corr = bad_data_double; 
      else           cnt_info[i_interp][i].sp_corr = ((count*fo_sum) - (f_sum*o_sum))/den;
      
      // Compute Kendall Tau Rank correlation coefficient:
      // For each pair of ranked data points (fi, oi), compare it to all other pairs
      // of ranked data points (fj, oj) where j > i.  If the relative ordering of the
      // ranks of the f's is the same as the relative ordering of the ranks of the o's,
      // count the comparison as concordant.  If the previous is not the case, count
      // the comparison as discordant.  If there is a tie between the o's, count the
      // comparison as extra_f.  A tie between the f's counts as an extra_o.  If there
      // is a tie in both the f's and o's, don't count the comparison as anything.
      concordant = discordant = extra_f = extra_o = 0;
      for(j=0; j<n_fcst_rank; j++) {
         for(k=j+1; k<n_fcst_rank; k++) {

            if(j==k) continue;

            // Check for agreement in the relative ordering of ranks
            if(      (fcst_rank[j] > fcst_rank[k] && obs_rank[j] > obs_rank[k]) ||
                     (fcst_rank[j] < fcst_rank[k] && obs_rank[j] < obs_rank[k]) ) concordant++;
            // Check for disagreement in the relative ordering of ranks
            else if( (fcst_rank[j] > fcst_rank[k] && obs_rank[j] < obs_rank[k]) ||
                     (fcst_rank[j] < fcst_rank[k] && obs_rank[j] > obs_rank[k]) ) discordant++;
            // Check for ties in the forecast rank
            else if(  fcst_rank[j] == fcst_rank[k] && obs_rank[j] != obs_rank[k]) extra_o++;
            // Check for ties in the observation rank
            else if(  fcst_rank[j] != fcst_rank[k] && obs_rank[j] == obs_rank[k]) extra_f++;
         }
      }
      den = sqrt((double) concordant+discordant+extra_f)*
            sqrt((double) concordant+discordant+extra_o);
      if(den == 0.0) cnt_info[i_interp][i].kt_corr = bad_data_double; 
      else           cnt_info[i_interp][i].kt_corr = (concordant - discordant)/den;
      
      // Compute the percntiles of the error
      sort(err_array, count);
      cnt_info[i_interp][i].e10 = percentile(err_array, count, 0.10);
      cnt_info[i_interp][i].e25 = percentile(err_array, count, 0.25);
      cnt_info[i_interp][i].e50 = percentile(err_array, count, 0.50);
      cnt_info[i_interp][i].e75 = percentile(err_array, count, 0.75);
      cnt_info[i_interp][i].e90 = percentile(err_array, count, 0.90);

      // Compute mean error and standard deviation of the mean error
      cnt_info[i_interp][i].me.v = err_sum/count;
      cnt_info[i_interp][i].estdev.v = compute_stdev(err_sum, err_sq_sum, count);
            
      // Compute mean absolute error
      cnt_info[i_interp][i].mae = abs_err_sum/count;
            
      // Compute the mean squared error
      cnt_info[i_interp][i].mse = err_sq_sum/count;
      
      // Compute the bias corrected mean squared error (decomposition of MSE)
      f = cnt_info[i_interp][i].fbar.v;
      o = cnt_info[i_interp][i].obar.v;
      cnt_info[i_interp][i].bcmse = cnt_info[i_interp][i].mse - (f-o)*(f-o);      
      
      // Compute root mean squared error
      cnt_info[i_interp][i].rmse = sqrt(err_sq_sum/count);      

      // Compute confidence intervals for each alpha value
      // In computing the confidence intervals, the spatial correlation between
      // adjacent points is ignored, and certain assumptions of normality are
      // made.  The user must decide if the computation method is appropriate
      // for the chosen field.
      for(j=0; j<cnt_info[i_interp][i].n_alpha; j++) {
            
         if(verbosity > 1) {
            cout << "[" << field_name
                 << ", " << mask_names[i]
                 << ", " << conf.ci_alpha(j).dval() << "]\n" << flush;
         }

         // Store the alpha value in the CNTInfo object
         cnt_info[i_interp][i].alpha[j] = conf.ci_alpha(j).dval();
         
         // Check for the degenerate case
         if(count <= 1) {
            cnt_info[i_interp][i].fbar.v_cl[j]    = cnt_info[i_interp][i].fbar.v_cu[j]    = bad_data_double;
            cnt_info[i_interp][i].fstdev.v_cl[j]  = cnt_info[i_interp][i].fstdev.v_cu[j]  = bad_data_double;
            cnt_info[i_interp][i].obar.v_cl[j]    = cnt_info[i_interp][i].obar.v_cu[j]    = bad_data_double;
            cnt_info[i_interp][i].ostdev.v_cl[j]  = cnt_info[i_interp][i].ostdev.v_cu[j]  = bad_data_double;
            cnt_info[i_interp][i].pr_corr.v_cl[j] = cnt_info[i_interp][i].pr_corr.v_cu[j] = bad_data_double;
            cnt_info[i_interp][i].me.v_cl[j]      = cnt_info[i_interp][i].me.v_cu[j]      = bad_data_double;
            cnt_info[i_interp][i].estdev.v_cl[j]  = cnt_info[i_interp][i].estdev.v_cu[j]  = bad_data_double;
      
            continue;
         }

         // Compute the critical values for the Normal or Student's-T distribution
         // based on the sample size
         if(count >= large_sample_threshold) {
            cv_normal_l = normal_cdf_inv(cnt_info[i_interp][i].alpha[j]/2.0, 0.0, 1.0);
            cv_normal_u = normal_cdf_inv(1.0 - (cnt_info[i_interp][i].alpha[j]/2.0), 0.0, 1.0);
         }
         // If the number of samples is less than the large sample threshold, 
         // use the T-distribution
         else {
            cv_normal_l = students_t_cdf_inv(cnt_info[i_interp][i].alpha[j]/2.0, count - 1);
            cv_normal_u = students_t_cdf_inv(1.0 - (cnt_info[i_interp][i].alpha[j]/2.0), count - 1);
         }

         // Compute the critical values for the Chi Squared distribution
         cv_chi2_l = chi2_cdf_inv(cnt_info[i_interp][i].alpha[j]/2.0, count - 1);
         cv_chi2_u = chi2_cdf_inv(1.0 - (cnt_info[i_interp][i].alpha[j]/2.0), count - 1);

         // Compute confidence interval for forecast mean
         cnt_info[i_interp][i].fbar.v_cl[j] = cnt_info[i_interp][i].fbar.v + 
                                    cv_normal_l*cnt_info[i_interp][i].fstdev.v/sqrt((double) count);
         cnt_info[i_interp][i].fbar.v_cu[j] = cnt_info[i_interp][i].fbar.v +
                                    cv_normal_u*cnt_info[i_interp][i].fstdev.v/sqrt((double) count);

         // Compute confidence interval for forecast standard deviation,
         // assuming normality of the forecast values
         cnt_info[i_interp][i].fstdev.v_cl[j] = sqrt((count - 1)*cnt_info[i_interp][i].fstdev.v*cnt_info[i_interp][i].fstdev.v/cv_chi2_u);
         cnt_info[i_interp][i].fstdev.v_cu[j] = sqrt((count - 1)*cnt_info[i_interp][i].fstdev.v*cnt_info[i_interp][i].fstdev.v/cv_chi2_l);

         // Compute confidence interval for observation mean
         cnt_info[i_interp][i].obar.v_cl[j] = cnt_info[i_interp][i].obar.v + 
                                    cv_normal_l*cnt_info[i_interp][i].ostdev.v/sqrt((double) count);
         cnt_info[i_interp][i].obar.v_cu[j] = cnt_info[i_interp][i].obar.v + 
                                    cv_normal_u*cnt_info[i_interp][i].ostdev.v/sqrt((double) count);

         // Compute confidence interval for observation standard deviation
         // assuming normality of the observation values
         cnt_info[i_interp][i].ostdev.v_cl[j] = sqrt((count - 1)*cnt_info[i_interp][i].ostdev.v*cnt_info[i_interp][i].ostdev.v/cv_chi2_u);
         cnt_info[i_interp][i].ostdev.v_cu[j] = sqrt((count - 1)*cnt_info[i_interp][i].ostdev.v*cnt_info[i_interp][i].ostdev.v/cv_chi2_l);

         // Compute confidence interval for the Pearson correlation coefficient
         if(cnt_info[i_interp][i].pr_corr.v == bad_data_double || count <= 3) {
            cnt_info[i_interp][i].pr_corr.v_cl[j] = bad_data_double;
            cnt_info[i_interp][i].pr_corr.v_cu[j] = bad_data_double;
         }
         else {
            v = 0.5*log((1 + cnt_info[i_interp][i].pr_corr.v)/(1 - cnt_info[i_interp][i].pr_corr.v));
            cl = v + cv_normal_l/sqrt((double) (count - 3));
            cu = v + cv_normal_u/sqrt((double) (count - 3));
            cnt_info[i_interp][i].pr_corr.v_cl[j] = (pow(e, 2*cl) - 1)/(pow(e, 2*cl) + 1);
            cnt_info[i_interp][i].pr_corr.v_cu[j] = (pow(e, 2*cu) - 1)/(pow(e, 2*cu) + 1);
         }

         // Compute confidence interval for mean error
         cnt_info[i_interp][i].me.v_cl[j] = cnt_info[i_interp][i].me.v + 
                                  cv_normal_l*cnt_info[i_interp][i].estdev.v/sqrt((double) count);
         cnt_info[i_interp][i].me.v_cu[j] = cnt_info[i_interp][i].me.v + 
                                  cv_normal_u*cnt_info[i_interp][i].estdev.v/sqrt((double) count);

         // Compute confidence interval for the error standard deviation
         cnt_info[i_interp][i].estdev.v_cl[j] = sqrt((count - 1)*cnt_info[i_interp][i].estdev.v*cnt_info[i_interp][i].estdev.v/cv_chi2_u);
         cnt_info[i_interp][i].estdev.v_cu[j] = sqrt((count - 1)*cnt_info[i_interp][i].estdev.v*cnt_info[i_interp][i].estdev.v/cv_chi2_l);
      } // end for j
   } // end for i

   if(err_array) { delete [] err_array; err_array = (double *) 0; }
   if(fcst_rank) { delete [] fcst_rank; fcst_rank = (double *) 0; }
   if(obs_rank)  { delete [] obs_rank;  obs_rank  = (double *) 0; }
   
   return;
}

////////////////////////////////////////////////////////////////////////

void write_fcst_netcdf(const WrfData &fcst_wd, const WrfData &obs_wd, 
                       int lev, int mthd, int wdth) {
   int i, n, x, y, status, mon, day, yr, hr, min, sec;
   float *fcst_data = (float *) 0;
   float *diff_data = (float *) 0;
   double fval, oval;
   unixtime ut;
   NcVar *fcst_var, *diff_var;
   char diff_var_name[max_str_len], fcst_var_name[max_str_len];
   char time_str[max_str_len];
   
   // Allocate memory for the forecast, observation, and difference fields
   fcst_data = new float [grid.nx()*grid.ny()];
   diff_data = new float [grid.nx()*grid.ny()];

   // Compute the difference field for each of the masking regions
   for(i=0; i<n_masks; i++) {

      // If the neighborhood width is greater than 1, name the variables
      // using the smoothing method
      if(wdth > 1) {
         sprintf(fcst_var_name, "FCST_%s_%s_%s_%s_%i", 
                 get_grib_code_abbr(gc_info[lev].code, conf.ncep_defaults().ival(), status),
                 gc_info[lev].lvl_str, mask_names[i], mthd_str[mthd], wdth*wdth);
         sprintf(diff_var_name, "DIFF_%s_%s_%s_%s_%i", 
                 get_grib_code_abbr(gc_info[lev].code, conf.ncep_defaults().ival(), status),
                 gc_info[lev].lvl_str, mask_names[i], mthd_str[mthd], wdth*wdth);
      }
      // If the neighborhood width is equal to 1, no smoothing was done
      else {
         sprintf(fcst_var_name, "FCST_%s_%s_%s", 
                 get_grib_code_abbr(gc_info[lev].code, conf.ncep_defaults().ival(), status),
                 gc_info[lev].lvl_str, mask_names[i]);
         sprintf(diff_var_name, "DIFF_%s_%s_%s", 
                 get_grib_code_abbr(gc_info[lev].code, conf.ncep_defaults().ival(), status),
                 gc_info[lev].lvl_str, mask_names[i]);
      }

      // Define the forecast and difference variables
      fcst_var = nc_out->add_var(fcst_var_name, ncFloat, lat_dim, lon_dim);
      diff_var = nc_out->add_var(diff_var_name, ncFloat, lat_dim, lon_dim);
   
      // Add variable attributes for the forecast field
      fcst_var->add_att("name", field_name);
      fcst_var->add_att("long_name", get_grib_code_name(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      fcst_var->add_att("units", get_grib_code_unit(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      fcst_var->add_att("level", level_name);

      ut = fcst_wd.get_valid_time();
      unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
      sprintf(time_str, "%.4i-%.2i-%.2i %.2i:%.2i:%.2i", yr, mon, day, hr, min, sec);
      fcst_var->add_att("valid", time_str);

      fcst_var->add_att("masking_region", mask_names[i]);
      fcst_var->add_att("smoothing_method", mthd_str[mthd]);
      fcst_var->add_att("smoothing_neighborhood", wdth*wdth);
      fcst_var->add_att("_FillValue", bad_data_float);
            
      // Add variable attributes for the difference field
      diff_var->add_att("name", field_name);
      diff_var->add_att("long_name", get_grib_code_name(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      diff_var->add_att("units", get_grib_code_unit(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      diff_var->add_att("level", level_name);

      ut = fcst_wd.get_valid_time();
      unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
      sprintf(time_str, "%.4i-%.2i-%.2i %.2i:%.2i:%.2i", yr, mon, day, hr, min, sec);
      diff_var->add_att("fcst_valid", time_str);
      
      ut = obs_wd.get_valid_time();
      unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
      sprintf(time_str, "%.4i-%.2i-%.2i %.2i:%.2i:%.2i", yr, mon, day, hr, min, sec);
      diff_var->add_att("obs_valid", time_str);      

      diff_var->add_att("masking_region", mask_names[i]);
      diff_var->add_att("smoothing_method", mthd_str[mthd]);
      diff_var->add_att("smoothing_neighborhood", wdth*wdth);
      diff_var->add_att("_FillValue", bad_data_float);

      // Store the forecast, and difference values
      for(x=0; x<grid.nx(); x++) {
         for(y=0; y<grid.ny(); y++) {

            n = two_to_one(grid.nx(), grid.ny(), x, y);

            // Check whether the mask is on or off for this point
            if(!mask_wd[i].s_is_on(x, y)) {
               fcst_data[n] = diff_data[n] = bad_data_float;
               continue;
            }

            // Retrieve the forecast and observation values
            fval = fcst_wd.get_xy_double(x, y);
            oval = obs_wd.get_xy_double(x, y);

            // Set the forecast data
            if(fval == bad_data_double) fcst_data[n] = bad_data_float;
            else                        fcst_data[n] = fval;

            // Set the difference data
            if(fval == bad_data_double ||
               oval == bad_data_double) diff_data[n] = bad_data_float;
            else                        diff_data[n] = fval - oval;
         } // end for y
      } // end for x
   
      // Write out the forecast field
      if(!fcst_var->put(&fcst_data[0], grid.ny(), grid.nx())) {
         cerr << "\n\nERROR: write_fcst_netcdf() -> "
              << "error with the fcst_var->put for grib code " << field_name
              << " and masking region " << mask_names[i] << "\n\n" << flush;
         exit(1);
      }
      
      // Write out the difference field
      if(!diff_var->put(&diff_data[0], grid.ny(), grid.nx())) {
         cerr << "\n\nERROR: write_fcst_netcdf() -> "
              << "error with the diff_var->put for grib code " << field_name
              << " and masking region " << mask_names[i] << "\n\n" << flush;
         exit(1);
      }
   
   } // end for i
   
   // Deallocate and clean up
   if(fcst_data) { delete [] fcst_data; fcst_data = (float *) 0; }
   if(diff_data) { delete [] diff_data; diff_data = (float *) 0; }

   return;
}

////////////////////////////////////////////////////////////////////////

void write_obs_netcdf(const WrfData &obs_wd, int lev) {
   int i, n, x, y, status, mon, day, yr, hr, min, sec;
   float *obs_data = (float *) 0;
   double oval;
   unixtime ut;
   NcVar *obs_var;
   char time_str[max_str_len], obs_var_name[max_str_len];
   
   // Allocate memory for the observation field
   obs_data  = new float [grid.nx()*grid.ny()];

   // Compute the difference field for each of the masking regions
   for(i=0; i<n_masks; i++) {

      sprintf(obs_var_name, "OBS_%s_%s_%s", 
              get_grib_code_abbr(gc_info[lev].code, conf.ncep_defaults().ival(), status),
              gc_info[lev].lvl_str, mask_names[i]);

      // Define the observation variable
      obs_var  = nc_out->add_var(obs_var_name, ncFloat, lat_dim, lon_dim);
   
      // Add variable attributes for the observation field
      obs_var->add_att("name", field_name);
      obs_var->add_att("long_name", get_grib_code_name(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      obs_var->add_att("units", get_grib_code_unit(gc_info[lev].code, conf.ncep_defaults().ival(), status));
      obs_var->add_att("level", level_name);
      
      ut = obs_wd.get_valid_time();
      unix_to_mdyhms(ut, mon, day, yr, hr, min, sec);
      sprintf(time_str, "%.4i-%.2i-%.2i %.2i:%.2i:%.2i", yr, mon, day, hr, min, sec);
      obs_var->add_att("valid", time_str);      

      obs_var->add_att("masking_region", mask_names[i]);
      obs_var->add_att("_FillValue", bad_data_float);

      // Store the forecast, observation, and difference values
      for(x=0; x<grid.nx(); x++) {
         for(y=0; y<grid.ny(); y++) {

            n = two_to_one(grid.nx(), grid.ny(), x, y);

            // Check whether the mask is on or off for this point
            if(!mask_wd[i].s_is_on(x, y)) {
               obs_data[n] = bad_data_float;
               continue;
            }

            // Retrieve the observation value
            oval = obs_wd.get_xy_double(x, y);

            // Set the observation data
            if(oval == bad_data_double) obs_data[n]  = bad_data_float;
            else                        obs_data[n]  = oval;
         } // end for y
      } // end for x
      
      // Write out the observation field
      if(!obs_var->put(&obs_data[0], grid.ny(), grid.nx())) {
         cerr << "\n\nERROR: write_obs_netcdf() -> "
              << "error with the obs_var->put for grib code " << field_name
              << " and masking region " << mask_names[i] << "\n\n" << flush;
         exit(1);
      }
   } // end for i
   
   // Deallocate and clean up
   if(obs_data)  { delete [] obs_data;  obs_data  = (float *) 0; }

   return;
}

////////////////////////////////////////////////////////////////////////

int get_interp_mthd(const char *str) {
   int i, n;

   n = -1;
   
   for(i=0; i<max_mthd; i++) {
      if(strcmp(mthd_str[i], str) == 0) {
         n = i;
         break;
      }
   }
   
   if(n == -1) {
      cerr << "\n\nERROR: get_interp_mthd() -> "
           << "Can't find matching interpolation method for string \""
           << str << "\"!\n\n" << flush;
      exit(1);
   }
   
   return(n);
}

////////////////////////////////////////////////////////////////////////

WrfData smooth_field(WrfData *wd, int mthd, int wdth) {
   WrfData smooth_wd;
   double v;
   int x, y, x_ll, y_ll;
   
   // Initialize the smoothed field to the raw field
   smooth_wd = *wd;
   
   // Check for a width value of 1 for which no smoothing should
   // be performed
   if(wdth <= 1) return(smooth_wd);

   // Otherwise, apply smoothing to each grid point
   for(x=0; x<wd->get_nx(); x++) {
      for(y=0; y<wd->get_ny(); y++) {
      
         // The neighborhood width must be odd, find the lower-left 
         // corner of the neighborhood
         x_ll = x - (wdth - 1)/2;
         y_ll = y - (wdth - 1)/2;

         // Compute the smoothed value based on the interpolation method
         switch(mthd) {

            case i_min:     // Minimum
               v = interp_min(wd, x_ll, y_ll, wdth, 
                              conf.interp_threshold().dval());
               break;
      
            case i_max:     // Maximum
               v = interp_max(wd, x_ll, y_ll, wdth,
                              conf.interp_threshold().dval());
               break;
      
            case i_median:  // Median
               v = interp_median(wd, x_ll, y_ll, wdth,
                                 conf.interp_threshold().dval());
               break;
      
            case i_uw_mean: // Unweighted Mean
               v = interp_uw_mean(wd, x_ll, y_ll, wdth,
                                  conf.interp_threshold().dval());
               break;

            // Distance-weighted mean is omitted here since it is not
            // an option for grid_stat

            // Least-squares fit is omitted here since it is not
            // an option for grid_stat
 
            default:
               cerr << "\n\nERROR: smooth_field() -> "
                    << "unexpected interpolation method encountered: "
                    << mthd << "\n\n" << flush;
               exit(1);
               break;
         }

         // Store the smoothed value
         smooth_wd.put_xy_double(v, x, y);
      } // end for y
   } // end for x
   
   return(smooth_wd);
}

////////////////////////////////////////////////////////////////////////

void usage(int argc, char *argv[]) {
   
   cout << "\nUsage: " 
        << program_name << "\n"
        << "\tfcst_file\n"
        << "\tobs_file\n"
        << "\tconfig_file\n"
        << "\t[-outdir path]\n"
        << "\t[-v level]\n\n"

        << "\twhere\t\"fcst_file\" is a forecast file in either Grib "
        << "or netCDF format (output of pcp_combine) containing the "
        << "field(s) to be verified (required).\n"

        << "\t\t\"obs_file\" is an observation file in either Grib "
        << "or netCDF format (output of pcp_combine) containing the "
        << "verifying field(s) (required).\n"

        << "\t\t\"config_file\" is a GridStatConfig file containing "
        << "the desired configuration settings (required).\n"

        << "\t\t\"-outdir path\" overrides the default output directory (" 
        << out_dir << ") (optional).\n"

        << "\t\t\"-v level\" overrides the default level of logging ("
        << verbosity << ") (optional).\n"
 
        << "\n\tNOTE: The forecast and observation fields must be "
        << "on the same grid.\n\n" << flush;

   return;
}

////////////////////////////////////////////////////////////////////////