// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
// ** 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:   point_stat.cc
//
//   Description:
//      Based on user specified parameters, this tool compares a
//      a gridded forecast field to the point observation output of
//      the PB2NC or ASCII2NC tool.  It computes many verification
//      scores and statistics, including confidence intervals, to
//      summarize the comparison.
//
//   Mod#   Date      Name           Description
//   ----   ----      ----           -----------
//   000    04/18/07  Halley Gotway  New
//   001    12/20/07  Halley Gotway  Allow verification for level 0
//                    for verifying PRMSL
//   002    01/24/08  Halley Gotway  In compute_cnt, print a warning
//                    message when bad data values are encountered.
//   003    01/24/08  Halley Gotway  Add support for writing the matched
//                    pair data to the MPR file.
//   004    02/04/08  Halley Gotway  Modify to read new format of the
//                    intermediate NetCDF point-observation format.
//   005    02/11/08  Halley Gotway  Remove BIAS from the CNT line
//          since it's the same as the ME.
//   006    02/12/08  Halley Gotway  Fix bug in writing COP line to
//                    write out OY and ON rather than FY and FN.
//   007    02/12/08  Halley Gotway  Enable Point-Stat to read the
//                    NetCDF output of PCP-Combine.
//   008    02/25/08  Halley Gotway  Write the VSDB lines using the
//                    routines in the vx_met_util library.
//   009    07/01/08  Halley Gotway   Add the rank_corr_flag to the
//                    config file to disable computing rank correlations.
//   010    07/18/08  Halley Gotway   Fix bug in the write_mpr routine.
//                    Replace i_fho with i_mpr.
//
////////////////////////////////////////////////////////////////////////

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 "point_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"

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

static const char * program_name = "point_stat";

// Pointer to the random number generator to be used
static gsl_rng *rng_ptr = (gsl_rng *) 0;

static const int max_num_recs = 300;

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

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

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

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

// Input climo Grib or NetCDF files
static char        climo_file[max_str_len];
static GribFile    climo_gb_file;
static NcFile     *climo_nc_file;

// Input observation NetCDF file
static NcFile     *obs_in;
static StringArray obs_file;

static GribRecord fcst_r, climo_r;

// 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_sal1l2    = 8;
static const int i_vl1l2     = 9;
static const int i_val1l2    = 10;
static const int i_mpr       = 11;

static const int hdr_arr_len = 3; // Observation header length
static const int obs_arr_len = 5; // Observation values length

// Static variables for configuration file entries
static Grid grid;

// Flag to indicate the presence of a climatologly file
static int climo_flag = 0;

// 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 climo_fmt_flag = 0;

// Flag to indicate whether any vector winds will be
// computed
static int vector_flag = 0;

// Strings to be output in the VSDB file
static char abbr_str[max_str_len];
static char level_name[max_str_len];

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

// Objects to hold PairData information for each Grib Code
// Sized as gc_pd[n_gc]
static GCPairData *gc_pd;

// Objects to hold threshold information for each Grib Code
static ThreshInfo *gc_ti;

// Numbers of verifying message types, masks, and interpolation
// techniques
static int n_msg_typ     = 0;
static int n_mask        = 0;
static int n_mask_area   = 0;
static int n_mask_sid    = 0;
static int n_interp      = 0;
static int n_interp_wdth = 0;
static int n_interp_mthd = 0;

// Multiple verification masking regions to be applied to the forecast
// and observation fields
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;

// Array of WrfData objects for storing the forecast and climotological
// fields to be verified.
// Sized as fcst_wd[n_gc][n_fcst] and climo_wd[n_gc][n_climo]
static WrfData **fcst_wd, **climo_wd;

// Output file streams for creating VSDB and text files
static ofstream *vsdb_out, *cts_out, *cnt_out;
static ofstream *fho_out, *ctc_out, *ctp_out, *cfp_out, *cop_out;
static ofstream *sl1l2_out, *sal1l2_out, *vl1l2_out, *val1l2_out;
static ofstream *mpr_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];
static char sl1l2_file[max_str_len], sal1l2_file[max_str_len];
static char vl1l2_file[max_str_len], val1l2_file[max_str_len];
static char mpr_file[max_str_len];

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

static void process_command_line(int, char **);
static void process_fcst_climo_files();
static void process_config();
static void process_grid(WrfData &);
static void process_masks();
static void parse_sid_mask(const char *, StringArray &);
static void setup_gc_pd_objects();
static void setup_txt_files(unixtime, int);
static void process_grib_codes();
static void process_obs_file(int);
static void process_scores();
static void clean_up();

static void compute_cnt_cts(CNTInfo &, CTSInfo *&, int, PairData *);
static void compute_sl1l2(SL1L2Info &, int, PairData *);
static void compute_vl1l2(VL1L2Info &, int, PairData *, int, PairData *);

static void write_cnt(CNTInfo &, int, PairData *);
static void write_cts(CTSInfo &, int, PairData *);
static void write_sl1l2(SL1L2Info &, int, PairData *);
static void write_vl1l2(VL1L2Info &, int, PairData *, int, PairData *);
static void write_mpr(int, PairData *);
static void close_txt_files();

static void clear_sl1l2(SL1L2Info &);
static void clear_vl1l2(VL1L2Info &);

static void usage(int, char **);

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

int main(int argc, char *argv[]) {
   WrfData wd;
   int i;

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

   // Process the forecast and climo files
   process_fcst_climo_files();

   // Process the configuration file
   process_config();

   // Process the grid to be used for this data
   process_grid(wd);

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

   // Setup the PairData objects
   setup_gc_pd_objects();

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

   // Process each of the grib codes to be verified
   process_grib_codes();

   // Process each observation netCDF file
   for(i=0; i<obs_file.n_elements(); i++) {
      process_obs_file(i);
   }

   // Compute the scores and write them out
   process_scores();

   // 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/point_stat", MET_BASE);

   if(argc < 4) {
      usage(argc, argv);
      exit(1);
   }

   strcpy(climo_file, "none");

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

   // Parse command line arguments
   for(i=0; i<argc; i++) {

      if(strcmp(argv[i], "-ncfile") == 0) {
         obs_file.add(argv[i+1]);
         i++;
      }
      else if(strcmp(argv[i], "-climo") == 0) {
         strcpy(climo_file, argv[i+1]);
         climo_flag = 1;
         i++;
      }
      else 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);

   // List the input files
   if(verbosity > 0) {
      cout << "Forecast File:\t" << fcst_file << "\n"
           << "Configuration File:\t" << config_file << "\n"
           << "Climatology File:\t" << climo_file << "\n" << flush;

      for(i=0; i<obs_file.n_elements(); i++) {
         cout << "Observation File:\t" << obs_file[i] << "\n" << flush;
      }
   }

   return;
}

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

void process_fcst_climo_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_climo_files() -> "
              << "can't open grib file: "
              << fcst_file << "\n\n" << flush;
         exit(1);
      }
   }

   // Open the climo file if specified
   if(climo_flag) {

      // Attempt to open the climo file as a NetCDF file.  If it
      // fails, attempt to read it as Grib.
      climo_nc_file = new NcFile(climo_file);

      if(!climo_nc_file || !(climo_nc_file->is_valid())) {
         climo_nc_file->close();
         delete climo_nc_file;
         climo_nc_file = (NcFile *) 0;
         climo_fmt_flag = 1;
      }
      else {
         climo_fmt_flag = 0;
      }

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

         if(!(climo_gb_file.open(climo_file)) ) {
            cerr << "\n\nERROR: process_fcst_climo_files() -> "
                 << "can't open grib file: "
                 << climo_file << "\n\n" << flush;
            exit(1);
         }
      }
   }

   return;
}

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

void process_config() {
   int i, n;
   GCInfo gc;

   //
   // Conf: model
   //

   if(strlen(conf.model().sval()) == 0 ||
      check_reg_exp(non_ws_reg_exp, conf.model().sval()) != true) {

      cerr << "\n\nERROR: process_config() -> "
           << "The model name must be set to a non-empty "
           << "string value.\n\n" << flush;
      exit(1);
   }

   //
   // 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_pd = new GCPairData [n_gc];

   for(i=0; i<n_gc; i++) {
      parse_gcinfo(conf.vx_grib_code(i).sval(), gc);

      // No support for WDIR
      if(gc.code == wdir_grib_code) {
         cerr << "\n\nERROR: process_config() -> "
              << "the wind direction field may not be verified using point_stat.\n\n"
              << flush;
         exit(1);
      }

      // Only support vertical levels at 0m, 2m, or 10m
      if(gc.lvl_type == VertLevel &&
         gc.lvl_1 != 0 &&
         gc.lvl_1 != 2 &&
         gc.lvl_1 != 10) {

         cerr << "\n\nERROR: process_config() -> "
              << "Only vertical level values of 2 and 10 meters "
              << "are currently supported \""
              << conf.vx_grib_code(i).sval() << "\"\n\n" << flush;
         exit(1);
      }
      gc_pd[i].set_gc(gc);
   }

   // Check the specified grib codes and turn on the vector_flag
   // when UGRD is followed VGRD at the same level
   for(i=0; i<n_gc; i++) {

      if(i+1 < n_gc
      && gc_pd[i].gc.code   == ugrd_grib_code
      && gc_pd[i+1].gc.code == vgrd_grib_code
      && strcmp(gc_pd[i].gc.lvl_str, gc_pd[i+1].gc.lvl_str) == 0) {

         gc_pd[i].gc.vector_flag   = 1;
         gc_pd[i+1].gc.vector_flag = 1;
         vector_flag = 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: message_type
   //

   // Check that at least one PrepBufr message type is provided
   if((n_msg_typ = conf.n_message_type_elements()) <= 0) {
      cerr << "\n\nERROR: process_config() -> "
           << "At least one PrepBufr message type must be provided.\n\n"
           << flush;
      exit(1);
   }

   // Check that each PrepBufr message type provided is valid
   for(i=0; i<n_msg_typ; i++) {

      if(strstr(all_msg_typ_str, conf.message_type(i).sval()) == NULL) {
         cerr << "\n\nERROR: process_config() -> "
              << "Invalid message type string provided ("
              << conf.message_type(i).sval() << ").\n\n" << flush;
         exit(1);
      }
   }

   //
   // 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: boot_interval
   //

   // Check that boot_interval is set to 0 or 1
   if(conf.boot_interval().ival() != boot_bca_flag &&
      conf.boot_interval().ival() != boot_perc_flag) {
      cerr << "\n\nERROR: process_config() -> "
           << "The boot_interval parameter must be set to "
           << boot_bca_flag << " or "
           << boot_perc_flag << "!\n\n" << flush;
      exit(1);
   }

   //
   // Conf: boot_rep_prop
   //

   // Check that boot_rep_prop is set between 0 and 1
   if(conf.boot_rep_prop().dval() <= 0.0 ||
      conf.boot_rep_prop().dval() > 1.0) {
      cerr << "\n\nERROR: process_config() -> "
           << "The boot_rep_prop parameter must be set between "
           << "0 and 1!\n\n" << flush;
      exit(1);
   }

   //
   // Conf: n_boot_rep
   //

   // Check that n_boot_rep is set > 0
   if(conf.n_boot_rep().dval() < 0.0) {
      cerr << "\n\nERROR: process_config() -> "
           << "The number of bootstrap resamples (n_boot_rep) "
           << "must be set to a value >= 0.\n\n" << flush;
      exit(1);
   }

   //
   // Conf: boot_rng and boot_seed
   //

   // Set the random number generator and seed value to be used when
   // computing bootstrap confidence intervals
   rng_set(rng_ptr, conf.boot_rng().sval(), conf.boot_seed().sval());

   //
   // Conf: interp_method
   //

   // Check that at least one interpolation method is provided
   if((n_interp_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);
   }

   // Get the interpolation method to make sure it's valid
   for(i=0; i<n_interp_mthd; i++) {
      n = get_interp_mthd(conf.interp_method(i).sval());
   }

   //
   // Conf: interp_width
   //

   // Check that at least one interpolation width is provided
   if((n_interp_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_interp_wdth; i++) {

      if(conf.interp_width(i).ival() < 1) {
         cerr << "\n\nERROR: process_config() -> "
              << "The interpolation width values must be set "
              << "greater than or equal to 1 ("
              << conf.interp_width(i).ival() << ").\n\n" << flush;
         exit(1);
      }

      if(conf.interp_width(i).ival() == 1) n_interp += 1;
      if(conf.interp_width(i).ival() >  1) n_interp += n_interp_mthd;
   }

   //
   // 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);
   }

   // Store the interpolation threshold value
   for(i=0; i<n_gc; i++) {
      gc_pd[i].set_interp_threshold(conf.interp_threshold().dval());
   }

   //
   // Conf: output_flag
   //

   // Check that at least one output VSDB type is requested
   for(i=i_fho, n=0; i<i_val1l2; i++) n += (conf.output_flag(i).ival() > 0);
   if(n == 0) {
      cerr << "\n\nERROR: process_config() -> "
           << "At least one output VSDB type must be requested.\n\n"
           << flush;
      exit(1);
   }

   //
   // Conf: rank_corr_flag
   //
   if(conf.rank_corr_flag().ival() != 0 &&
      conf.rank_corr_flag().ival() != 1) {
      cerr << "\n\nERROR: process_config() -> "
           << "The rank_corr_flag (" << conf.rank_corr_flag().ival()
           << ") must be set to 0 or 1.\n\n"
           << flush;
      exit(1);
   }

   return;
}

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

void process_grid(WrfData &wd) {
   Grid fcst_grid;

   // Process the forecast file in NetCDF format (output of pcp_combine)
   if(fcst_fmt_flag == 0) {

      // Get the grib code abbreviation string
      get_grib_code_abbr(gc_pd[0].gc.code, conf.ncep_defaults().ival(),
                         abbr_str);

      read_pcp_combine_netcdf(fcst_nc_file, abbr_str, level_name,
                              wd, fcst_grid, verbosity,
                              pc_data, lc_data, st_data);
   }
   // Process the forecast file in Grib format
   else {

      // Read the first grib record from the forecast file to extract
      // the grid information
      read_grib_record(fcst_gb_file, fcst_r, 0, wd, fcst_grid, verbosity,
                       pc_data, lc_data, st_data);
   }

   // Store the grid
   grid = fcst_grid;

   return;
}

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

void process_masks() {
   int i, n_mask_grid, n_mask_poly;
   StringArray mask_sid;

   // Get the number of masking areas
   n_mask_grid = conf.n_mask_grids_elements();
   n_mask_poly = conf.n_mask_polys_elements();
   n_mask_area = n_mask_grid + n_mask_poly;

   // Get the number of masking station ID's
   parse_sid_mask(conf.mask_stations().sval(), mask_sid);
   n_mask_sid = mask_sid.n_elements();

   // Save the total number masks as a sum of the masking areas and
   // masking points
   n_mask = n_mask_area + n_mask_sid;

   // Check that at least one verification masking region is provided
   if(n_mask == 0) {
      cerr << "\n\nERROR: process_masks() -> "
           << "At least one grid, polyline, or station ID "
           << "masking region must be provided.\n\n" << flush;
      exit(1);
   }

   // Allocate space to store the masking information
   mask_wd    = new WrfData [n_mask_area];
   mask_names = new char   *[n_mask];

   // Parse out the masking grid areas
   for(i=0; i<n_mask_grid; i++)
      parse_grid_mask(conf.mask_grids(i).sval(), grid,
                      mask_wd[i],
                      mask_names[i]);

   // Parse out the masking poly areas
   for(i=0; i<n_mask_poly; i++)
      parse_poly_mask(conf.mask_polys(i).sval(), grid,
                      mask_wd[i+n_mask_grid],
                      mask_names[i+n_mask_grid]);

   // Store the masking station ID points
   for(i=0; i<n_mask_sid; i++) {
      mask_names[i+n_mask_grid+n_mask_poly] = new char [strlen(mask_sid[i])+1];
      strcpy(mask_names[i+n_mask_grid+n_mask_poly], mask_sid[i]);
   }

   return;
}

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

void parse_sid_mask(const char *mask_sid_file, StringArray &mask_sid) {
   ifstream in;
   char mask_sid_file_tmp[max_str_len], sid_str[max_str_len];

   // Replace any instances of MET_BASE with it's expanded value
   replace_string("MET_BASE", MET_BASE, mask_sid_file, mask_sid_file_tmp);

   // Check for an empty length string
   if(strlen(mask_sid_file_tmp) == 0) return;

   // Open the mask station id file specified
   in.open(mask_sid_file_tmp);

   if(!in) {
      cerr << "\n\nERROR: parse_sid_mask() -> "
           << "Can't open the mask station ID file specified ("
           << mask_sid_file_tmp << ").\n\n" << flush;
      exit(1);
   }

   // Read in and store the masking station ID names
   while(in >> sid_str) mask_sid.add(sid_str);

   // Close the input file
   in.close();

   return;
}

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

void setup_gc_pd_objects() {
   int i, j, k, n;

   // PairData is stored in the gc_pd objects in the following order:
   // [n_msg_typ][n_mask][n_interp]
   for(i=0; i<n_gc; i++) {

      // Set up the dimensions for the gc_pd object
      gc_pd[i].set_pd_size(n_msg_typ, n_mask, n_interp);

      // Add the verifying message type to the gc_pd objects
      for(j=0; j<n_msg_typ; j++)
         gc_pd[i].set_msg_typ(j, conf.message_type(j).sval());

      // Add the masking information to the gc_pd objects
      for(j=0; j<n_mask; j++) {

         // If this is a masking area
         if(j<n_mask_area)
            gc_pd[i].set_mask_wd(j, mask_names[j], &mask_wd[j]);
         // Otherwise this is a masking StationID
         else
            gc_pd[i].set_mask_wd(j, mask_names[j], (WrfData *) 0);
      }

      // Add the interpolation methods to the gc_pd objects
      for(j=0, n=0; j<n_interp_wdth; j++) {

         // Check for width equal to 1
         if(conf.interp_width(j).ival() == 1) {
            gc_pd[i].set_interp(n, i_uw_mean, 1);
            n += 1;
         }
         else if(conf.interp_width(j).ival() > 1) {

            for(k=0; k<n_interp_mthd; k++) {
               gc_pd[i].set_interp(n, conf.interp_method(k).sval(),
                                   conf.interp_width(j).ival());
               n += 1;
            } // end for k
         }
      } // end for j
   } // end for i

   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], tmp_str[max_str_len];

   // Initialize output file streams
   vsdb_out  = fho_out    = cts_out   = cnt_out    = (ofstream *) 0;
   ctc_out   = ctp_out    = cfp_out   = cop_out    = (ofstream *) 0;
   sl1l2_out = sal1l2_out = vl1l2_out = val1l2_out = (ofstream *) 0;
   mpr_out                                         = (ofstream *) 0;

   // Construct 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(tmp_str, "%s/point_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(vsdb_file,   "%s.vsdb",       tmp_str);
   sprintf(fho_file,    "%s_fho.txt",    tmp_str);
   sprintf(ctc_file,    "%s_ctc.txt",    tmp_str);
   sprintf(ctp_file,    "%s_ctp.txt",    tmp_str);
   sprintf(cfp_file,    "%s_cfp.txt",    tmp_str);
   sprintf(cop_file,    "%s_cop.txt",    tmp_str);
   sprintf(cts_file,    "%s_cts.txt",    tmp_str);
   sprintf(cnt_file,    "%s_cnt.txt",    tmp_str);
   sprintf(sl1l2_file,  "%s_sl1l2.txt",  tmp_str);
   sprintf(sal1l2_file, "%s_sal1l2.txt", tmp_str);
   sprintf(vl1l2_file,  "%s_vl1l2.txt",  tmp_str);
   sprintf(val1l2_file, "%s_val1l2.txt", tmp_str);
   sprintf(mpr_file,    "%s_mpr.txt",    tmp_str);

   // 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:\t"
           << 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:\t"
              << fho_file << "\n" << flush;
      }

      get_line_names(fho_columns, n_fho_columns, 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:\t"
              << ctc_file << "\n" << flush;
      }

      get_line_names(ctc_columns, n_ctc_columns, 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:\t"
              << ctp_file << "\n" << flush;
      }

      get_line_names(ctp_columns, n_ctp_columns, 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:\t"
              << cfp_file << "\n" << flush;
      }

      get_line_names(cfp_columns, n_cfp_columns, 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:\t"
              << cop_file << "\n" << flush;
      }

      get_line_names(cop_columns, n_cop_columns, 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:\t"
              << cts_file << "\n" << flush;
      }

      get_line_names(cts_columns, n_cts_columns, 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:\t"
              << cnt_file << "\n" << flush;
      }

      get_line_names(cnt_columns, n_cnt_columns, 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:\t"
              << sl1l2_file << "\n" << flush;
      }

      get_line_names(sl1l2_columns, n_sl1l2_columns, header);
      *sl1l2_out << header << "\n" << flush;
   }

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

      // Create and open the SAL1L2 output file stream
      sal1l2_out = new ofstream;
      sal1l2_out->open(sal1l2_file);

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

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

      get_line_names(sal1l2_columns, n_sal1l2_columns, header);
      *sal1l2_out << header << "\n" << flush;
   }

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

      // Create and open the VL1L2 output file stream
      vl1l2_out = new ofstream;
      vl1l2_out->open(vl1l2_file);

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

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

      get_line_names(vl1l2_columns, n_vl1l2_columns, header);
      *vl1l2_out << header << "\n" << flush;
   }

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

      // Create and open the VAL1L2 output file stream
      val1l2_out = new ofstream;
      val1l2_out->open(val1l2_file);

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

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

      get_line_names(val1l2_columns, n_val1l2_columns, header);
      *val1l2_out << header << "\n" << flush;
   }

   // Create the matched pairs text file
   if(conf.output_flag(i_mpr).ival() >= 2) {

      // Create and open the matched pairs output file stream
      mpr_out = new ofstream;
      mpr_out->open(mpr_file);

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

      // List the matched pairs file as it is being created
      if(verbosity > 1) {
         cout << "Creating Output MPR file:\t"
              << mpr_file << "\n" << flush;
      }

      get_line_names(mpr_columns, n_mpr_columns, header);
      *mpr_out << header << "\n" << flush;
   }

   return;
}

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

void process_grib_codes() {
   int i, j, status;
   int n_fcst_rec,  fcst_rec[max_num_recs], fcst_lvl[max_num_recs];
   int n_climo_rec, climo_rec[max_num_recs], climo_lvl[max_num_recs];
   Grid fcst_grid, climo_grid;
   char tmp_str[max_str_len];

   // Allocate space to store the forecast and climo fields for each
   // grib code.
   fcst_wd  = new WrfData * [n_gc];
   climo_wd = new WrfData * [n_gc];

   // Loop through each of the grib codes to be verified and extract
   // the forecast and climatological fields from the grib files needed
   // for verification
   for(i=0; i<n_gc; i++) {

      // Initialize the WrfData pointers
      fcst_wd[i]  = (WrfData *) 0;
      climo_wd[i] = (WrfData *) 0;

      // Store the grib code abbreviation in a string
      get_grib_code_abbr(gc_pd[i].gc.code, conf.ncep_defaults().ival(),
                         abbr_str);

      if(verbosity > 1 && gc_pd[i].gc.lvl_str == (char *) 0) {
         cout << "\n----------------------------------------\n\n"
              << "Reading records for Grib Field " << abbr_str
              << "\n" << flush;
      }
      else if(verbosity > 1) {
         cout << "\n----------------------------------------\n\n"
              << "Reading records for Grib Field " << abbr_str
              << " at " << gc_pd[i].gc.lvl_str << "\n" << flush;
      }

      // Process the forecast file in NetCDF format (output of pcp_combine)
      if(fcst_fmt_flag == 0) {
         n_fcst_rec = 1;
      }
      // Grib format
      // Compute a list of grib records for this grib code which
      // fall within the requested range of levels from the forecast
      // and climatological data, including one level above and one
      // level below the range
      else {
         n_fcst_rec = find_grib_record_levels(fcst_gb_file, gc_pd[i].gc.code,
                                              gc_pd[i].gc.lvl_1, gc_pd[i].gc.lvl_2,
                                              fcst_rec, fcst_lvl);

         if(n_fcst_rec == 0) {
            cout << "***WARNING***: process_grib_codes() -> "
                 << "no records matching grib code " << gc_pd[i].gc.code
                 << " with level indicator of "
                 << gc_pd[i].gc.lvl_str << " not found in grib file: "
                 << fcst_file << "\n" << flush;
            continue;
         }
      }

      // Process the climo file if specified
      if(climo_flag) {

         // NetCDF format
         if(climo_fmt_flag == 0) {
            n_climo_rec = 1;
         }
         // Grib format
         else {
            n_climo_rec = find_grib_record_levels(climo_gb_file, gc_pd[i].gc.code,
                                                  gc_pd[i].gc.lvl_1, gc_pd[i].gc.lvl_2,
                                                  climo_rec, climo_lvl);

            if(n_climo_rec == 0) {
               cout << "***WARNING***: process_grib_codes() -> "
                    << "no records matching grib code " << gc_pd[i].gc.code
                    << " with level indicator of "
                    << gc_pd[i].gc.lvl_str << " not found in grib file: "
                    << climo_file << "\n" << flush;
               continue;
            }
         }
      }
      // No climo file specified
      else {
         n_climo_rec = 0;
      }

      // Dump out the number of levels found
      if(verbosity > 1) {
         cout << "For Grib Field " << abbr_str << " at "
           << gc_pd[i].gc.lvl_str << " found " << n_fcst_rec
           << " forecast levels and " << n_climo_rec
           << " climatology levels\n" << flush;
      }

      // Allocate space to store the forecast and climo fields
      if(n_fcst_rec > 0) fcst_wd[i]  = new WrfData [n_fcst_rec];
      if(n_fcst_rec > 0) climo_wd[i] = new WrfData [n_climo_rec];

      // Set the number of pointers to the raw forecast and climo fields
      gc_pd[i].set_n_fcst(n_fcst_rec);
      gc_pd[i].set_n_climo(n_climo_rec);

      // Read in the forecast fields for this grib code
      for(j=0; j<n_fcst_rec; j++) {

         // NetCDF format
         if(fcst_fmt_flag == 0) {
            read_pcp_combine_netcdf(fcst_nc_file, abbr_str, tmp_str,
                                    fcst_wd[i][j], fcst_grid, verbosity,
                                    pc_data, lc_data, st_data);
         }
         // Grib format
         else {
            status = get_grib_record(fcst_gb_file, fcst_r,
                                     gc_pd[i].gc.code, fcst_lvl[j], fcst_lvl[j],
                                     fcst_wd[i][j], fcst_grid, verbosity,
                                     pc_data, lc_data, st_data);

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

         // Check that the grid dimensions have not changed
         if(fcst_grid.nx() != grid.nx() ||
            fcst_grid.ny() != grid.ny()) {
            cerr << "\n\nERROR: process_grib_codes() -> "
                 << "The forecast grid dimensions "
                 << "have changed (" << fcst_grid.nx() << ", "
                 << fcst_grid.ny() << ") != (" << grid.nx()
                 << ", " << grid.ny() << ") for grib code "
                 << gc_pd[i].gc.code << ".\n\n" << flush;
            exit(1);
         }

         // Store information for the raw forecast fields
         gc_pd[i].set_fcst_prs_lvl(j, fcst_lvl[j]);
         gc_pd[i].set_fcst_wd_ptr(j, &fcst_wd[i][j]);

      } // end for j

      // Read in the climatology fields for this grib code
      for(j=0; j<n_climo_rec; j++) {

         // NetCDF format
         if(climo_fmt_flag == 0) {
            read_pcp_combine_netcdf(climo_nc_file, abbr_str, tmp_str,
                                    climo_wd[i][j], climo_grid, verbosity,
                                    pc_data, lc_data, st_data);
         }
         // Grib format
         else {
            status = get_grib_record(climo_gb_file, climo_r,
                                     gc_pd[i].gc.code, climo_lvl[j], climo_lvl[j],
                                     climo_wd[i][j], climo_grid, verbosity,
                                     pc_data, lc_data, st_data);

            if(status != 0) {
               cout << "***WARNING***: process_grib_codes() -> "
                    << "grib code " << gc_pd[i].gc.code
                    << " with accumulation/level indicator of "
                    << climo_lvl[j] << " not found in grib file: "
                    << climo_file << "\n" << flush;
               continue;
            }
         }

         // Check that the grid dimensions have not changed
         if(climo_grid.nx() != grid.nx() ||
            climo_grid.ny() != grid.ny()) {
            cerr << "\n\nERROR: process_grib_codes() -> "
                 << "The climatology grid dimensions "
                 << "have changed (" << climo_grid.nx() << ", "
                 << climo_grid.ny() << ") != (" << grid.nx()
                 << ", " << grid.ny() << ") for grib code "
                 << gc_pd[i].gc.code << ".\n\n" << flush;
            exit(1);
         }

         // Store information for the raw climo fields
         gc_pd[i].set_climo_prs_lvl(j, climo_lvl[j]);
         gc_pd[i].set_climo_wd_ptr(j, &climo_wd[i][j]);

      } // end for j
   } // end for i

   if(verbosity > 1) {
      cout << "\n----------------------------------------\n\n" << flush;
   }

   return;
}

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

void process_obs_file(int i_nc) {
   int i, j, k, l, i_obs;
   float obs_arr[obs_arr_len], hdr_arr[hdr_arr_len];
   float prev_obs_arr[obs_arr_len];
   char hdr_typ_str[max_str_len], hdr_sid_str[max_str_len];
   PairData *pd_ptr;

   // Open the observation file as a NetCDF file.
   // The observation file must be in NetCDF format as the
   // output of the PB2NC or ASCII2NC tool.
   obs_in = new NcFile(obs_file[i_nc]);

   if(!obs_in || !(obs_in->is_valid())) {
      obs_in->close();
      delete obs_in;
      obs_in = (NcFile *) 0;

      cerr << "\n\nERROR: process_obs_file() -> "
           << "can't open observation netCDF file: "
           << obs_file[i_nc] << "\n\n" << flush;
      exit(1);
   }

   // Define dimensions
   NcDim *strl_dim    = (NcDim *) 0; // Maximum string length
   NcDim *obs_dim     = (NcDim *) 0; // Number of observations
   NcDim *hdr_dim     = (NcDim *) 0; // Number of PrepBufr messages

   // Define variables
   NcVar *obs_arr_var = (NcVar *) 0;
   NcVar *hdr_typ_var = (NcVar *) 0;
   NcVar *hdr_sid_var = (NcVar *) 0;
   NcVar *hdr_vld_var = (NcVar *) 0;
   NcVar *hdr_arr_var = (NcVar *) 0;

   // Read the dimensions
   strl_dim = obs_in->get_dim("mxstr");
   obs_dim  = obs_in->get_dim("nobs");
   hdr_dim  = obs_in->get_dim("nhdr");

   if(!strl_dim || !strl_dim->is_valid() ||
      !obs_dim  || !obs_dim->is_valid()  ||
      !hdr_dim  || !hdr_dim->is_valid()) {
      cerr << "\n\nERROR: process_obs_file() -> "
           << "can't read \"mxstr\", \"nobs\" or \"nmsg\" "
           << "dimensions from netCDF file: "
           << obs_file[i_nc] << "\n\n" << flush;
      exit(1);
   }

   // Read the variables
   obs_arr_var = obs_in->get_var("obs_arr");
   hdr_typ_var = obs_in->get_var("hdr_typ");
   hdr_sid_var = obs_in->get_var("hdr_sid");
   hdr_vld_var = obs_in->get_var("hdr_vld");
   hdr_arr_var = obs_in->get_var("hdr_arr");

   if(!obs_arr_var || !obs_arr_var->is_valid() ||
      !hdr_typ_var || !hdr_typ_var->is_valid() ||
      !hdr_sid_var || !hdr_sid_var->is_valid() ||
      !hdr_vld_var || !hdr_vld_var->is_valid() ||
      !hdr_arr_var || !hdr_arr_var->is_valid()) {
      cerr << "\n\nERROR: process_obs_file() -> "
           << "can't read \"obs_arr\", \"hdr_typ\", \"hdr_sid\", "
           << "\"hdr_vld\", or \"hdr_arr\" variables from netCDF file: "
           << obs_file[i_nc] << "\n\n" << flush;
      exit(1);
   }

   if(verbosity > 1) {
      cout << "Processing " << obs_dim->size()
           << " observations from " << hdr_dim->size()
           << " PrepBufr messages...\n" << flush;
   }

   // Process each observation in the file
   for(i_obs=0; i_obs<obs_dim->size(); i_obs++) {

      // Read the current observation message
      if(!obs_arr_var->set_cur((long) i_obs)
      || !obs_arr_var->get(obs_arr, 1, obs_arr_len)) {
         cerr << "\n\nERROR: process_obs_file() -> "
              << "can't read the record for observation "
              << "number " << i_obs << "\n\n" << flush;
         exit(1);
      }

      // If the current observation is UGRD, save it as the
      // previous.  If vector winds are to be computed, UGRD
      // must be followed by VGRD
      if(vector_flag && obs_arr[1] == ugrd_grib_code) {
         for(j=0; j<4; j++) prev_obs_arr[j] = obs_arr[j];
      }

      // If the current observation is VGRD and vector
      // winds are to be computed.  Make sure that the
      // previous observation was UGRD with the same header
      // and at the same vertical level.
      if(vector_flag && obs_arr[1] == vgrd_grib_code) {

         if(obs_arr[0] != prev_obs_arr[0] ||
            obs_arr[2] != prev_obs_arr[2] ||
            obs_arr[3] != prev_obs_arr[3]) {
            cerr << "\n\nERROR: process_obs_file() -> "
                 << "when computing vector winds, each UGRD observation "
                 << "must be followed by a VGRD observation with the same "
                 << "header and at the same level\n\n" << flush;
            exit(1);
         }
      }

      // Read the corresponding header array for this observation
      if(!hdr_arr_var->set_cur((long) (obs_arr[0]-1)) ||
         !hdr_arr_var->get(hdr_arr, 1, hdr_arr_len)) {
         cerr << "\n\nERROR: process_obs_file() -> "
              << "can't read the header array record for header "
              << "number " << obs_arr[0]-1 << "\n\n" << flush;
         exit(1);
      }

      // Read the corresponding header type for this observation
      if(!hdr_typ_var->set_cur((long) (obs_arr[0]-1)) ||
         !hdr_typ_var->get(hdr_typ_str, 1, strl_dim->size())) {
         cerr << "\n\nERROR: process_obs_file() -> "
              << "can't read the message type record for header "
              << "number " << obs_arr[0]-1 << "\n\n" << flush;
         exit(1);
      }

      // Read the corresponding header Station ID for this observation
      if(!hdr_sid_var->set_cur((long) (obs_arr[0]-1)) ||
         !hdr_sid_var->get(hdr_sid_str, 1, strl_dim->size())) {
         cerr << "\n\nERROR: process_obs_file() -> "
              << "can't read the station ID record for header "
              << "number " << obs_arr[0]-1 << "\n\n" << flush;
         exit(1);
      }

      // Check each gc_pd object to see if this observation should be added
      for(j=0; j<n_gc; j++) {

         // Attempt to add the observation to the gc_pd object
         gc_pd[j].add_obs(hdr_arr, hdr_typ_str, hdr_sid_str, obs_arr, grid);
      }
   } // end for i_obs

   // Dump out the number of pairs found
   if(verbosity > 1) {
      cout << "\n";
      for(i=0; i<n_gc; i++) {
         for(j=0; j<gc_pd[i].n_msg_typ; j++) {
            for(k=0; k<gc_pd[i].n_mask; k++) {
               for(l=0; l<gc_pd[i].n_interp; l++) {

                  get_grib_code_abbr(gc_pd[i].gc.code,
                                     conf.ncep_defaults().ival(),
                                     abbr_str);

                  pd_ptr = &gc_pd[i].pd[j][k][l];

                  cout << "For Grib Field " << abbr_str
                       << " at " << gc_pd[i].gc.lvl_str
                       << " for observation type " << pd_ptr->msg_typ
                       << ", over region " << pd_ptr->mask_name
                       << ", for interpolation method " << mthd_str[pd_ptr->interp_mthd]
                       << "(" << pd_ptr->interp_wdth << "), found "
                       << pd_ptr->n_pair << " pairs.\n" << flush;
               } // end for l
            } // end for k
         }// end for j
      } // end for i
      cout << "\n";
   } // end if

   // Deallocate and clean up
   obs_in->close();
   delete obs_in;
   obs_in = (NcFile *) 0;

   return;
}

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

void process_scores() {
   int i, j, k, l, m;
   PairData *pd_ptr;
   CTSInfo  *cts_info;
   CNTInfo   cnt_info;
   SL1L2Info sl1l2_info;
   VL1L2Info vl1l2_info;

   // Loop through the Grib Code Pair Data objects, compute the CTS,
   // CNT, SL1L2, and SAL1L2 quantities, and write the output.
   for(i=0; i<n_gc; i++) {
      for(j=0; j<gc_pd[i].n_msg_typ; j++) {
         for(k=0; k<gc_pd[i].n_mask; k++) {
            for(l=0; l<gc_pd[i].n_interp; l++) {

               pd_ptr = &gc_pd[i].pd[j][k][l];

               // Continue if the number of points is <= 1
               if(pd_ptr->n_pair == 0) continue;

               // Compute CNT and CTS scores for each threshold specified
               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() ||
                  conf.output_flag(i_cnt).ival()) {

                  // Allocate memory for CTSInfo objects
                  cts_info = new CTSInfo [gc_ti[i].n_thresh];

                  cnt_info.clear();
                  compute_cnt_cts(cnt_info, cts_info, i, pd_ptr);

                  // Write out the CNT statistics
                  if(conf.output_flag(i_cnt).ival()) {
                     write_cnt(cnt_info, i, pd_ptr);
                  } // end write CNT

                  // Write out the CTS statistics
                  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()) {

                     for(m=0; m<gc_ti[i].n_thresh; m++) {
                        write_cts(cts_info[m], i, pd_ptr);
                     }
                  } // end write CTS

                  // Deallocate memory for CTSInfo objects
                  if(cts_info) { delete [] cts_info; cts_info = (CTSInfo *) 0; }

               } // end compute CNT and CTS

               // Compute SL1L2 and SAL1L2 scores as long as the vector_flag is not set
               if((conf.output_flag(i_sl1l2).ival() || conf.output_flag(i_sal1l2).ival())
               && gc_pd[i].gc.vector_flag == 0) {
                  clear_sl1l2(sl1l2_info);
                  compute_sl1l2(sl1l2_info, i, pd_ptr);
                  write_sl1l2(sl1l2_info, i, pd_ptr);
               }

               // Compute VL1L2 and VAL1L2 scores for the VGRD grib code
               // preceeded by the UGRD grib code
               if((conf.output_flag(i_vl1l2).ival() || conf.output_flag(i_val1l2).ival())
               && gc_pd[i].gc.vector_flag == 1
               && i > 0
               && gc_pd[i].gc.code == vgrd_grib_code
               && gc_pd[i-1].gc.code == ugrd_grib_code) {
                  clear_vl1l2(vl1l2_info);
                  compute_vl1l2(vl1l2_info, i-1, &gc_pd[i-1].pd[j][k][l],
                                i, &gc_pd[i].pd[j][k][l]);
                  write_vl1l2(vl1l2_info, i-1, &gc_pd[i-1].pd[j][k][l],
                              i, &gc_pd[i].pd[j][k][l]);
               }

               // Write out the MPR lines
               if(conf.output_flag(i_mpr).ival()) {
                  write_mpr(i, pd_ptr);
               }
            } // end for l
         } // end for k
      } // end for j

      if(verbosity > 1) {
         cout << "\n----------------------------------------\n\n" << flush;
      }
   } // end for i

   return;

}

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

void clean_up() {
   int i;

   // Close the output text files that were open for writing
   close_txt_files();

   // Close the forecast file
   if(fcst_fmt_flag == 0) fcst_nc_file->close();
   else                   fcst_gb_file.close();

   // Close the climatology file
   if(     climo_flag && climo_fmt_flag == 0) climo_nc_file->close();
   else if(climo_flag && climo_fmt_flag == 1) climo_gb_file.close();

   // Deallocate memory
   if(gc_pd)            { delete [] gc_pd;         gc_pd = (GCPairData *) 0; }
   if(gc_ti)            { delete [] gc_ti;         gc_ti = (ThreshInfo *) 0; }
   if(mask_wd)          { delete [] mask_wd;       mask_wd = (WrfData *) 0; }

   for(i=0; i<n_mask; i++) {
      if(mask_names[i]) { delete [] mask_names[i]; mask_names[i] = (char *) 0; }
   }
   if(mask_names)       { delete [] mask_names;    mask_names = (char **) 0; }

   for(i=0; i<n_gc; i++) {
      if(fcst_wd[i])    { delete [] fcst_wd[i];    fcst_wd[i] = (WrfData *) 0; }
      if(climo_wd[i])   { delete [] climo_wd[i];   climo_wd[i] = (WrfData *) 0; }
   }
   if(fcst_wd)          { delete [] fcst_wd;       fcst_wd = (WrfData **) 0; }
   if(climo_wd)         { delete [] climo_wd;      climo_wd = (WrfData **) 0; }

   // Deallocate memory for the random number generator
   rng_free(rng_ptr);

   return;
}

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

void compute_cnt_cts(CNTInfo &cnt_info, CTSInfo *&cts_info,
                     int i_gc, PairData *pd_ptr) {
   int i, j, n_cts;
   NumArray f_na, o_na;

   //
   // Set up the CNTInfo alpha values
   //
   cnt_info.allocate_n_alpha(n_ci_alpha);
   for(i=0; i<n_ci_alpha; i++) {
      cnt_info.alpha[i] = conf.ci_alpha(i).dval();
   }

   //
   // Set up the CTSInfo thresholds and alpha values
   //
   n_cts = gc_ti[i_gc].n_thresh;
   for(i=0; i<n_cts; i++) {
      cts_info[i].cts_thresh = gc_ti[i_gc].thresh[i];
      cts_info[i].allocate_n_alpha(n_ci_alpha);

      for(j=0; j<n_ci_alpha; j++) {
         cts_info[i].alpha[j] = conf.ci_alpha(j).dval();
      }
   }

   if(verbosity > 1) {

      get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                         abbr_str);

      cout << "Computing CNT and CTS for Grib Field " << abbr_str
           << " at " << gc_pd[i_gc].gc.lvl_str
           << " for observation type " << pd_ptr->msg_typ
           << ", over region " << pd_ptr->mask_name
           << ", for interpolation method " << mthd_str[pd_ptr->interp_mthd]
           << "(" << pd_ptr->interp_wdth << ") using "
           << pd_ptr->n_pair << " pairs\n" << flush;
   }

   //
   // Keep track of the interpolation method and width
   //
   cnt_info.interp_mthd = pd_ptr->interp_mthd;
   cnt_info.interp_wdth = pd_ptr->interp_wdth;
   for(i=0; i<n_cts; i++) {
      cts_info[i].interp_mthd = pd_ptr->interp_mthd;
      cts_info[i].interp_wdth = pd_ptr->interp_wdth;
   }

   //
   // Store the pairs in NumArray objects
   //
   for(i=0; i<pd_ptr->n_pair; i++) {
      f_na.add(pd_ptr->fdata[i]);
      o_na.add(pd_ptr->odata[i]);
   }

   //
   // Compute the stats, normal confidence intervals, and bootstrap
   // bootstrap confidence intervals
   //
   if(conf.boot_interval().ival() == boot_bca_flag) {
      compute_stats_ci_bca(rng_ptr, f_na, o_na, gc_pd[i_gc].gc.code,
         conf.n_boot_rep().ival(),
         cnt_info, cts_info, n_cts, conf.rank_corr_flag().ival(),
         conf.tmp_dir().sval());
   }
   else {
      compute_stats_ci_perc(rng_ptr, f_na, o_na, gc_pd[i_gc].gc.code,
         conf.n_boot_rep().ival(), conf.boot_rep_prop().dval(),
         cnt_info, cts_info, n_cts, conf.rank_corr_flag().ival(),
         conf.tmp_dir().sval());
   }

   return;
}

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

void compute_sl1l2(SL1L2Info &s_info, int i_gc, PairData *pd_ptr) {
   int i;
   double f, o, c;
   double f_sum, o_sum, fo_sum, ff_sum, oo_sum;
   double fa_sum, oa_sum, foa_sum, ffa_sum, ooa_sum;

   if(verbosity > 1) {

      get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                         abbr_str);

      cout << "Computing SL1L2 for Grib Field " << abbr_str
           << " at " << gc_pd[i_gc].gc.lvl_str
           << " for observation type " << pd_ptr->msg_typ
           << ", over region " << pd_ptr->mask_name
           << ", for interpolation method " << mthd_str[pd_ptr->interp_mthd]
           << "(" << pd_ptr->interp_wdth << ") using "
           << pd_ptr->n_pair << " pairs\n" << flush;
   }

   // Keep track of the interpolation method and width
   s_info.interp_mthd = pd_ptr->interp_mthd;
   s_info.interp_wdth = pd_ptr->interp_wdth;

   // Initialize the counts and sums
   s_info.scount = s_info.sacount = 0;
   f_sum  = o_sum  = fo_sum  = ff_sum  = oo_sum  = 0;
   fa_sum = oa_sum = foa_sum = ffa_sum = ooa_sum = 0;

   // Loop through the pair data and compute sums
   for(i=0; i<pd_ptr->n_pair; i++) {

      f = pd_ptr->fdata[i];
      o = pd_ptr->odata[i];
      c = pd_ptr->cdata[i];

      // Skip bad data values in the forecast or observation fields
      if(f == bad_data_double || o == bad_data_double) continue;

      // SL1L2 sums
      f_sum  += f;
      o_sum  += o;
      fo_sum += f*o;
      ff_sum += f*f;
      oo_sum += o*o;

      s_info.scount++;

      // SAL1L2 sums
      if(c != bad_data_double) {

         fa_sum  += f-c;
         oa_sum  += o-c;
         foa_sum += (f-c)*(o-c);
         ffa_sum += (f-c)*(f-c);
         ooa_sum += (o-c)*(o-c);

         s_info.sacount++;
      }
   }

   if(s_info.scount == 0) {
      cerr << "\n\nERROR: compute_sl1l2() -> "
           << "count is zero!\n\n" << flush;
      exit(1);
   }

   // Compute the mean SL1L2 values
   s_info.fbar  = f_sum/s_info.scount;
   s_info.obar  = o_sum/s_info.scount;
   s_info.fobar = fo_sum/s_info.scount;
   s_info.ffbar = ff_sum/s_info.scount;
   s_info.oobar = oo_sum/s_info.scount;

   // Compute the mean SAL1L2 values
   if(s_info.sacount > 0) {
      s_info.fabar  = fa_sum/s_info.sacount;
      s_info.oabar  = oa_sum/s_info.sacount;
      s_info.foabar = foa_sum/s_info.sacount;
      s_info.ffabar = ffa_sum/s_info.sacount;
      s_info.ooabar = ooa_sum/s_info.sacount;
   }
   else {
      s_info.fabar  = bad_data_double;
      s_info.oabar  = bad_data_double;
      s_info.foabar = bad_data_double;
      s_info.ffabar = bad_data_double;
      s_info.ooabar = bad_data_double;
   }

   return;
}

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

void compute_vl1l2(VL1L2Info &v_info,
                   int i_ugrd, PairData *ugrd_pd_ptr,
                   int i_vgrd, PairData *vgrd_pd_ptr) {
   char abbr_str_2[max_str_len];
   int i;
   double uf, vf, uc, vc, uo, vo;
   double uf_sum, vf_sum, uo_sum, vo_sum, uvfo_sum, uvff_sum, uvoo_sum;
   double ufa_sum, vfa_sum, uoa_sum, voa_sum, uvfoa_sum, uvffa_sum, uvooa_sum;

   if(verbosity > 1) {

      get_grib_code_abbr(gc_pd[i_ugrd].gc.code, conf.ncep_defaults().ival(),
                         abbr_str);
      get_grib_code_abbr(gc_pd[i_vgrd].gc.code, conf.ncep_defaults().ival(),
                         abbr_str_2);

      cout << "Computing VL1L2 for Grib Field " << abbr_str << "/"
           << abbr_str_2
           << " at " << gc_pd[i_ugrd].gc.lvl_str
           << " for observation type " << ugrd_pd_ptr->msg_typ
           << ", over region " << ugrd_pd_ptr->mask_name
           << ", for interpolation method " << mthd_str[ugrd_pd_ptr->interp_mthd]
           << "(" << ugrd_pd_ptr->interp_wdth << ") using "
           << ugrd_pd_ptr->n_pair << " pairs\n" << flush;
   }

   // Check that the number of pairs are the same
   if(ugrd_pd_ptr->n_pair != vgrd_pd_ptr->n_pair) {
      cerr << "\n\nERROR: compute_vl1l2() -> "
           << "the number of ugrd pairs should equal the number "
           << "of vgrd pairs: " << ugrd_pd_ptr->n_pair
           << " != " << vgrd_pd_ptr->n_pair
           << "\n\n" << flush;
      exit(1);
   }

   // Keep track of the interpolation method and width
   v_info.interp_mthd = ugrd_pd_ptr->interp_mthd;
   v_info.interp_wdth = ugrd_pd_ptr->interp_wdth;

   // Initialize the counts and sums
   v_info.vcount = v_info.vacount = 0;
   uf_sum  = vf_sum  = uo_sum  = vo_sum  = uvfo_sum  = uvff_sum  = uvoo_sum  = 0.0;
   ufa_sum = vfa_sum = uoa_sum = voa_sum = uvfoa_sum = uvffa_sum = uvooa_sum = 0.0;

   // Loop through the pair data and compute sums
   for(i=0; i<ugrd_pd_ptr->n_pair; i++) {

      uf = ugrd_pd_ptr->fdata[i];
      vf = vgrd_pd_ptr->fdata[i];
      uc = ugrd_pd_ptr->cdata[i];
      vc = vgrd_pd_ptr->cdata[i];
      uo = ugrd_pd_ptr->odata[i];
      vo = vgrd_pd_ptr->odata[i];

      // Skip bad data values in the forecast or observation fields
      if(uf == bad_data_double || vf == bad_data_double ||
         uo == bad_data_double || vo == bad_data_double) continue;

      // VL1L2 sums
      uf_sum   += uf;
      vf_sum   += vf;
      uo_sum   += uo;
      vo_sum   += vo;
      uvfo_sum += uf*uo+vf*vo;
      uvff_sum += uf*uf+vf*vf;
      uvoo_sum += uo*uo+vo*vo;

      v_info.vcount++;

      // VAL1L2 sums
      if(uc != bad_data_double && vc != bad_data_double) {

         ufa_sum   += uf-uc;
         vfa_sum   += vf-vc;
         uoa_sum   += uo-uc;
         voa_sum   += vo-vc;
         uvfoa_sum += (uf-uc)*(uo-uc)+(vf-vc)*(vo-vc);
         uvffa_sum += (uf-uc)*(uf-uc)+(vf-vc)*(vf-vc);
         uvooa_sum += (uo-uc)*(uo-uc)+(vo-vc)*(vo-vc);

         v_info.vacount++;
      }
   }

   if(v_info.vcount == 0) {
      cerr << "\n\nERROR: compute_vl1l2() -> "
           << "count is zero!\n\n" << flush;
      exit(1);
   }

   // Compute the mean VL1L2 values
   v_info.ufbar   = uf_sum/v_info.vcount;
   v_info.vfbar   = vf_sum/v_info.vcount;
   v_info.uobar   = uo_sum/v_info.vcount;
   v_info.vobar   = vo_sum/v_info.vcount;
   v_info.uvfobar = uvfo_sum/v_info.vcount;
   v_info.uvffbar = uvff_sum/v_info.vcount;
   v_info.uvoobar = uvoo_sum/v_info.vcount;

   // Compute the mean VAL1L2 values
   if(v_info.vacount > 0) {
      v_info.ufabar   = ufa_sum/v_info.vacount;
      v_info.vfabar   = vfa_sum/v_info.vacount;
      v_info.uoabar   = uoa_sum/v_info.vacount;
      v_info.voabar   = voa_sum/v_info.vacount;
      v_info.uvfoabar = uvfoa_sum/v_info.vacount;
      v_info.uvffabar = uvffa_sum/v_info.vacount;
      v_info.uvooabar = uvooa_sum/v_info.vacount;
   }
   else {
      v_info.ufabar   = bad_data_double;
      v_info.vfabar   = bad_data_double;
      v_info.uoabar   = bad_data_double;
      v_info.voabar   = bad_data_double;
      v_info.uvfoabar = bad_data_double;
      v_info.uvffabar = bad_data_double;
      v_info.uvooabar = bad_data_double;
   }

   return;
}

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

void write_cnt(CNTInfo &cnt_info, int i_gc, PairData *pd_ptr) {
   char hdr_str[max_line_len], col_str[max_line_len];
   char type_str[max_str_len];
   char line[max_line_len];
   int i;

   //
   // Grib Field abbreviation
   //
   get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                      abbr_str);

   //
   // CNT line type
   //
   if(conf.output_flag(i_cnt).ival() != 0 &&
      cnt_info.n > 0) {

      for(i=0; i<n_ci_alpha; i++) {

         sprintf(type_str, "CNT/%.3f", cnt_info.alpha[i]);

         build_header_str(hdr_str,
            conf.model().sval(),
            gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
            gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
            pd_ptr->msg_typ,
            pd_ptr->mask_name,
            type_str,
            abbr_str,
            gc_pd[i_gc].gc.lvl_str);

         build_cnt_str(col_str, cnt_info, i);

         sprintf(line, "%s %s", hdr_str, col_str);

         *vsdb_out << line << "\n";

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

   return;
}

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

void write_cts(CTSInfo &cts_info, int i_gc, PairData *pd_ptr) {
   char hdr_str[max_line_len], col_str[max_line_len];
   char thresh_str[max_str_len], type_str[max_str_len];
   char line[max_line_len];
   int i;

   //
   // Type of thresholding performed
   //
   cts_info.cts_thresh.get_str(thresh_str);

   //
   // Grib Field abbreviation
   //
   get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                      abbr_str);

   //
   // FHO line type
   //
   if(conf.output_flag(i_fho).ival() != 0 &&
      cts_info.cts.n() > 0   ) {

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

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_fho_str(col_str, cts_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // CTC line type
   //
   if(conf.output_flag(i_ctc).ival() != 0 &&
      cts_info.cts.n() > 0) {

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

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_ctc_str(col_str, cts_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // CTP line type
   //
   if(conf.output_flag(i_ctp).ival() != 0 &&
      cts_info.cts.n() > 0) {

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

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_ctp_str(col_str, cts_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // CFP line type
   //
   if(conf.output_flag(i_cfp).ival() != 0 &&
      cts_info.cts.n() > 0) {

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

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_cfp_str(col_str, cts_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // COP line type
   //
   if(conf.output_flag(i_cop).ival() != 0 &&
      cts_info.cts.n() > 0) {

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

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_cop_str(col_str, cts_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // CTS line type
   //
   if(conf.output_flag(i_cts).ival() != 0 &&
      cts_info.cts.n() > 0) {

      for(i=0; i<cts_info.n_alpha; i++) {

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

         build_header_str(hdr_str,
            conf.model().sval(),
            gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
            gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
            pd_ptr->msg_typ,
            pd_ptr->mask_name,
            type_str,
            abbr_str,
            gc_pd[i_gc].gc.lvl_str);

         build_cts_str(col_str, cts_info, i);

         sprintf(line, "%s %s", hdr_str, col_str);

         *vsdb_out << line << "\n";

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

   return;
}

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

void write_sl1l2(SL1L2Info &s_info, int i_gc, PairData *pd_ptr) {
   char hdr_str[max_line_len], col_str[max_line_len];
   char type_str[max_str_len];
   char line[max_line_len];

   //
   // Grib Field abbreviation
   //
   get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                      abbr_str);

   //
   // SL1L2 line type
   //
   if(conf.output_flag(i_sl1l2).ival() != 0 &&
      s_info.scount > 0) {

      strcpy(type_str, "SL1L2");

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_sl1l2_str(col_str, s_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // SAL1L2 line type
   //
   if(conf.output_flag(i_sal1l2).ival() != 0 &&
      s_info.sacount > 0) {

      strcpy(type_str, "SAL1L2");

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      build_sal1l2_str(col_str, s_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   return;
}

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

void write_vl1l2(VL1L2Info &v_info,
                 int i_ugrd, PairData *ugrd_pd_ptr,
                 int i_vgrd, PairData *vgrd_pd_ptr) {
   char ugrd_str[max_str_len], vgrd_str[max_str_len];
   char hdr_str[max_line_len], col_str[max_line_len];
   char type_str[max_str_len];
   char line[max_line_len];

   //
   // Grib Field abbreviation
   //
   get_grib_code_abbr(gc_pd[i_ugrd].gc.code, conf.ncep_defaults().ival(),
                      ugrd_str);
   get_grib_code_abbr(gc_pd[i_vgrd].gc.code, conf.ncep_defaults().ival(),
                      vgrd_str);
   sprintf(abbr_str, "%s_%s", ugrd_str, vgrd_str);

   //
   // VL1L2 line type
   //
   if(conf.output_flag(i_vl1l2).ival() != 0 &&
      v_info.vcount > 0) {

      strcpy(type_str, "VL1L2");

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_ugrd].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_ugrd].fcst_wd_ptr[0]->get_valid_time(),
         ugrd_pd_ptr->msg_typ,
         ugrd_pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_ugrd].gc.lvl_str);

      build_vl1l2_str(col_str, v_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   //
   // VAL1L2 line type
   //
   if(conf.output_flag(i_val1l2).ival() != 0 &&
      v_info.vacount > 0) {

      strcpy(type_str, "VAL1L2");

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_ugrd].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_ugrd].fcst_wd_ptr[0]->get_valid_time(),
         ugrd_pd_ptr->msg_typ,
         ugrd_pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_ugrd].gc.lvl_str);

      build_val1l2_str(col_str, v_info);

      sprintf(line, "%s %s", hdr_str, col_str);

      *vsdb_out << line << "\n";

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

   return;
}

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

void write_mpr(int i_gc, PairData *pd_ptr) {
   char hdr_str[max_line_len], col_str[max_line_len];
   char type_str[max_str_len];
   char line[max_line_len];
   int i;

   //
   // Grib Field abbreviation
   //
   get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(),
                      abbr_str);

   //
   // MPR line type
   //
   if(conf.output_flag(i_mpr).ival() != 0 &&
      pd_ptr->n_pair > 0) {

      strcpy(type_str, "MPR");

      build_header_str(hdr_str,
         conf.model().sval(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_lead_time(),
         gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(),
         pd_ptr->msg_typ,
         pd_ptr->mask_name,
         type_str,
         abbr_str,
         gc_pd[i_gc].gc.lvl_str);

      for(i=0; i<pd_ptr->n_pair; i++) {

         build_mpr_str(col_str, pd_ptr, i);

         sprintf(line, "%s %s", hdr_str, col_str);

         *vsdb_out << line << "\n";

         if(conf.output_flag(i_mpr).ival() >= 2) *mpr_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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << sl1l2_file << "\n" << flush;
      }

      // Close the output SL1L2 file
      sl1l2_out->close();
      delete sl1l2_out;
      sl1l2_out = (ofstream *) 0;
   }

   // Close the output SAL1L2 txt file
   if(conf.output_flag(i_sal1l2).ival() >= 2 && sal1l2_out) {

      // List SAL1L2 file after it is finished
      if(verbosity > 0) {
         cout << "Output SAL1L2 file:\t"
              << sal1l2_file << "\n" << flush;
      }

      // Close the output SAL1L2 file
      sal1l2_out->close();
      delete sal1l2_out;
      sal1l2_out = (ofstream *) 0;
   }

   // Close the output VL1L2 txt file
   if(conf.output_flag(i_vl1l2).ival() >= 2 && vl1l2_out) {

      // List VL1L2 file after it is finished
      if(verbosity > 0) {
         cout << "Output VL1L2 file:\t"
              << vl1l2_file << "\n" << flush;
      }

      // Close the output VL1L2 file
      vl1l2_out->close();
      delete vl1l2_out;
      vl1l2_out = (ofstream *) 0;
   }

   // Close the output VAL1L2 txt file
   if(conf.output_flag(i_val1l2).ival() >= 2 && val1l2_out) {

      // List VAL1L2 file after it is finished
      if(verbosity > 0) {
         cout << "Output VAL1L2 file:\t"
              << val1l2_file << "\n" << flush;
      }

      // Close the output VAL1L2 file
      val1l2_out->close();
      delete val1l2_out;
      val1l2_out = (ofstream *) 0;
   }

   // Close the output matched pairs txt file
   if(conf.output_flag(i_mpr).ival() >= 2 && mpr_out) {

      // List the matched pairs file after it is finished
      if(verbosity > 0) {
         cout << "Output MPR file:\t"
              << mpr_file << "\n" << flush;
      }

      // Close the output Matched Pairs file
      mpr_out->close();
      delete mpr_out;
      mpr_out = (ofstream *) 0;
   }

   return;
}


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

void clear_sl1l2(SL1L2Info &s_info) {

   s_info.fbar    = 0.0;
   s_info.obar    = 0.0;
   s_info.fobar   = 0.0;
   s_info.ffbar   = 0.0;
   s_info.oobar   = 0.0;
   s_info.scount  = 0;

   s_info.fabar   = 0.0;
   s_info.oabar   = 0.0;
   s_info.foabar  = 0.0;
   s_info.ffabar  = 0.0;
   s_info.ooabar  = 0.0;
   s_info.sacount = 0;

   return;
}

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

void clear_vl1l2(VL1L2Info &v_info) {

   v_info.ufbar    = 0.0;
   v_info.vfbar    = 0.0;
   v_info.uobar    = 0.0;
   v_info.vobar    = 0.0;
   v_info.uvfobar  = 0.0;
   v_info.uvffbar  = 0.0;
   v_info.uvoobar  = 0.0;
   v_info.vcount   = 0;

   v_info.ufabar   = 0.0;
   v_info.vfabar   = 0.0;
   v_info.uoabar   = 0.0;
   v_info.voabar   = 0.0;
   v_info.uvfobar  = 0.0;
   v_info.uvffabar = 0.0;
   v_info.uvooabar = 0.0;
   v_info.vacount  = 0;

   return;
}

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

void usage(int argc, char *argv[]) {

   cout << "\nUsage: "
        << program_name << "\n"
        << "\tfcst_file\n"
        << "\tobs_file\n"
        << "\tconfig_file\n"
        << "\t[-climo climo_file]\n"
        << "\t[-ncfile netcdf_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 netCDF format "
        << "(output of PB2NC or ASCII2NC) containing the verifying "
        << "observation data points (required).\n"

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

        << "\t\t\"-climo climo_file\" is a climatological file in either "
        << "Grid or netCDF format (output of pcp_combine) on the same "
        << "grid as the forecast file to be\n"
        << "\t\tused when computing scalar and vector anomaly measures.  "
        << "If not provided, scalar and vector "
        << "anomaly values will not be computed (optional).\n"

        << "\t\t\"-ncfile netcdf_file\" may be used to specify "
        << "additional NetCDF point observation files to be used "
        << "(optional).\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"

        << flush;

   return;
}

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