// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
// ** 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 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    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 "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";

static const int max_num_vars = 50;
static const int max_num_recs = 300;
static const int var_str_len = 32;
static const int max_line_len = 1024;

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

// Variables for command line arguments
static point_stat_Conf conf;
static char fcst_file[max_str_len],   obs_file[max_str_len];
static char config_file[max_str_len], climo_file[max_str_len];
static char out_dir[max_str_len];
static int verbosity = 1;

// Variables for reading forecast and observation files
static GribFile   fcst_gb_file, climo_gb_file;
static GribRecord fcst_r, climo_r;
static NcFile    *obs_in;

// 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 hdr_arr_len = 7;     // PrepBufr message header length
static const int obs_arr_len = 11;    // Observation event array length

// Values of observation types in the PrepBufr file (P Q T Z U V)
static const int i_p_var  = 1;
static const int i_q_var  = 2;
static const int i_t_var  = 3;
static const int i_z_var  = 4;
static const int i_u_var  = 5;
static const int i_v_var  = 6;
// Derived observation types
static const int i_rh_var   = 7;
static const int i_cape_var = 8;
static const int i_cin_var  = 9;
static const int i_li_var   = 10;
static const int i_pbl_var  = 11;
static const int i_bndy_var = 12;

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

// Flag to indicate the presence of a climatologly file
static int climo_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 const char *met_version = "V01";
static char abbr_str[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 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 void process_command_line(int, char **);
static void process_fcst_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();
static void process_scores();
static void clean_up();

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

static void write_cts(CTSInfo &, int, PairData *);
static void write_cnt(CNTInfo &, int, PairData *);
static void write_sl1l2(SL1L2Info &, int, PairData *);
static void write_vl1l2(VL1L2Info &, int, PairData *, 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;
   
   // Process the command line arguments
   process_command_line(argc, argv);
   
   // Process the forecast and climo files
   process_fcst_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 the observation netCDF file
   process_obs_file();
   
   // 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]);
   strcpy(obs_file, argv[2]);
   strcpy(config_file, argv[3]);

   // Parse command line arguments
   for(i=0; i<argc; i++) {
   
      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: " << fcst_file << "\n" 
           << "Observation File: " << obs_file << "\n" 
           << "Configuration File: " << config_file << "\n"
           << "Climatology File: " << climo_file << "\n" << flush;
   }
   
   return;
}

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

void process_fcst_files() {

   // Open the forecast file as a Grib file.
   if( !(fcst_gb_file.open(fcst_file)) ) {
      cerr << "\n\nERROR: process_fcst_files() -> "
           << "can't open grib file: "
           << fcst_file << "\n\n" << flush;
      exit(1);
   }
   
   // If provided, open up the climatology file as a Grib file.
   if(climo_flag == 1) {

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

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

void process_config() {
   int i, n;
   GCInfo gc;
   
   //
   // Conf: vx_grib_code
   //
   
   // Parse out the grib codes to be verified
   if((n_gc = conf.n_vx_grib_code_elements()) <= 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 Pair Data objects
   gc_pd = new GCPairData [n_gc];
   
   // Parse out the Grib Code info and store it
   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: 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);
   }
   
   // Check that each interpolation method provided is valid
   for(i=0; i<n_interp_mthd; i++) {
   
      if(strstr(all_interp_mthd_str, conf.interp_method(i).sval()) == NULL) {
         cerr << "\n\nERROR: process_config() -> " 
              << "Invalid interpolation method string provided ("
              << conf.interp_method(i).sval() << ").\n\n" << flush;
         exit(1);
      }
   }

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

   return;
}

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

void process_grid(WrfData &wd) {
   Grid fcst_grid;
   
   // 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:
   // [mask][msg_typ][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];

   // 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;
   
   // 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/point_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/point_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/point_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/point_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/point_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/point_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/point_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/point_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/point_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);
   sprintf(sal1l2_file, "%s/point_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_sal1l2.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(vl1l2_file, "%s/point_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_vl1l2.txt",
           out_dir, l_hr, l_min, l_sec, yr, mon, day, hr, min, sec);
   sprintf(val1l2_file, "%s/point_stat_%.2i%.2i%.2iL_%.4i%.2i%.2i_%.2i%.2i%.2iV_val1l2.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:\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;
      }
      
      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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << 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:\t"
              << sl1l2_file << "\n" << flush;
      }
      
      sl1l2_hdr_str(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;
      }
      
      sal1l2_hdr_str(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;
      }
      
      vl1l2_hdr_str(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;
      }
      
      val1l2_hdr_str(header);
      *val1l2_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;

   // 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
      strcpy(abbr_str, get_grib_code_abbr(gc_pd[i].gc.code, 
             conf.ncep_defaults().ival(), status));
   
      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;         
      }
      
      // 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
      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(climo_flag) {
         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);
      }
      else {
         n_climo_rec = 0;
      }
      
      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;
      }

      if(climo_flag && 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;
      }
      
      // 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++) {

         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++) {

         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, j, k, l, i_obs, status;
   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 step.
   obs_in = new NcFile(obs_file);
   
   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 << "\n\n" << flush;
      exit(1);
   }

   // Define dimensions
   NcDim *strl_dim    = (NcDim *) 0; // Maximum string length
   NcDim *obs_dim     = (NcDim *) 0; // Number of observations
   NcDim *msg_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");
   msg_dim  = obs_in->get_dim("nmsg");
   
   if(!strl_dim || !strl_dim->is_valid() 
   || !obs_dim  || !obs_dim->is_valid()
   || !msg_dim  || !msg_dim->is_valid()) {
      cerr << "\n\nERROR: process_obs_file() -> "
           << "can't read \"mxstr\", \"nobs\" or \"nmsg\" "
           << "dimensions from netCDF file: "
           << obs_file << "\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 << "\n\n" << flush;
      exit(1);
   }
   
   if(verbosity > 1) {
      cout << "Processing " << obs_dim->size() 
           << " observations from " << msg_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[3] == ugrd_grib_code) {
         for(j=0; j<3; 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[3] == vgrd_grib_code) {
         
         if(obs_arr[0] != prev_obs_arr[0]
         || obs_arr[1] != prev_obs_arr[1]
         || obs_arr[2] != prev_obs_arr[2]) {
            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 data 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++) {

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

                  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
  
   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];

               // Compute 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())
                && pd_ptr->n_pair > 0) {

                  for(m=0; m<gc_ti[i].n_thresh; m++) {
                     cts_info.clear();
                     compute_cts(cts_info, i, pd_ptr, 
                                 gc_ti[i].thresh[m], gc_ti[i].thresh_ind[m]);
                     write_cts(cts_info, i, pd_ptr);
                  }
               }

               // Compute CNT scores
               if(conf.output_flag(i_cnt).ival()
               && pd_ptr->n_pair > 0) {
                  cnt_info.clear();
                  compute_cnt(cnt_info, i, pd_ptr);
                  write_cnt(cnt_info, i, pd_ptr);
               }

               // 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
               && pd_ptr->n_pair > 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
               && pd_ptr->n_pair > 0
               && 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]);
               }
            } // 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
   fcst_gb_file.close();
   
   // Close the climatology file
   if(climo_flag == 1) climo_gb_file.close();   

   // Close the netCDF observation file
   obs_in->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; }
   
   return;
}

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

void compute_cts(CTSInfo &cts_info, int i_gc, PairData *pd_ptr,
                 double t, int t_ind) {
   int i, j, status;
   double f, o;

   if(verbosity > 1) {
   
      strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(), status));
      
      cout << "Computing 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 << ") thresholding "
           << pd_ptr->n_pair << " pairs "
           << thresh_type_str[t_ind] << " " << t << "\n" << flush;
   }
   
   // Allocate enough space and store all of the alpha values to be used
   cts_info.allocate_n_alpha(n_ci_alpha);

   cts_info.cts_thresh = t;
   cts_info.cts_thresh_ind = t_ind;
   
   // Loop through the pair data and fill in the contingency table
   for(i=0; i<pd_ptr->n_pair; i++) {
   
      f = pd_ptr->fdata[i];
      o = pd_ptr->odata[i];
      
      switch(t_ind) {

         case thresh_lt:
            if(     f <  t && o <  t) cts_info.cts.inc_fy_oy();
            else if(f <  t && o >= t) cts_info.cts.inc_fy_on();
            else if(f >= t && o <  t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f >= t && o >= t 
            break;

         case thresh_le:
            if(     f <= t && o <= t) cts_info.cts.inc_fy_oy();
            else if(f <= t && o >  t) cts_info.cts.inc_fy_on();
            else if(f >  t && o <= t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f > t && o > t 
            break;

         case thresh_eq:
            if(     f == t && o == t) cts_info.cts.inc_fy_oy();
            else if(f == t && o != t) cts_info.cts.inc_fy_on();
            else if(f != t && o == t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f != t && o != t 
            break;

         case thresh_ne:
            if(     f != t && o != t) cts_info.cts.inc_fy_oy();
            else if(f != t && o == t) cts_info.cts.inc_fy_on();
            else if(f == t && o != t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f == t && o == t 
            break;

         case thresh_gt:
            if(     f >  t && o >  t) cts_info.cts.inc_fy_oy();
            else if(f >  t && o <= t) cts_info.cts.inc_fy_on();
            else if(f <= t && o >  t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f <= t && o <= t 
            break;

         case thresh_ge:
            if(     f >= t && o >= t) cts_info.cts.inc_fy_oy();
            else if(f >= t && o <  t) cts_info.cts.inc_fy_on();
            else if(f <  t && o >= t) cts_info.cts.inc_fn_oy();
            else                      cts_info.cts.inc_fn_on(); // f < t && o < t 
            break;

         default: 
            cerr << "\n\nERROR: compute_cts() -> "
                 << "unxpected threshold indicator value: "
                 << t_ind << "\n\n" << flush;
            exit(1);
            break;
      } // end switch
   } // end for i
   
   // Compute contingency table statistics
   cts_info.baser.v = cts_info.cts.oy_tp(status);
   cts_info.fmean.v = cts_info.cts.fy_tp(status);
   cts_info.acc.v   = cts_info.cts.accuracy(status);
   cts_info.fbias.v = cts_info.cts.fbias(status);
   cts_info.pody.v  = cts_info.cts.pod_yes(status);
   cts_info.podn.v  = cts_info.cts.pod_no(status);
   cts_info.pofd.v  = cts_info.cts.pofd(status);
   cts_info.far.v   = cts_info.cts.far(status);
   cts_info.csi.v   = cts_info.cts.csi(status);
   cts_info.gss.v   = cts_info.cts.gss(status);
   cts_info.hk.v    = cts_info.cts.hk(status);
   cts_info.hss.v   = cts_info.cts.hss(status);
   cts_info.odds.v  = cts_info.cts.odds(status);
      
   // For each choice of alpha, compute the corresponding 
   // confidence interval
   for(j=0; j<n_ci_alpha; j++) {

      // Store the alpha value in the CTSInfo object
      cts_info.alpha[j] = conf.ci_alpha(j).dval();

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

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

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

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

void compute_cnt(CNTInfo &cnt_info, int i_gc, PairData *pd_ptr) {
   int i, j, k, count, status;
   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;

   if(verbosity > 1) {
   
      strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(), status));
      
      cout << "Computing CNT 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;
   }
   
   // Allocate enough space and store all of the alpha values to be used
   cnt_info.allocate_n_alpha(n_ci_alpha);

   // Allocate space to store the differences for use in computing the 
   // percentiles and for storing the data ranks
   err_array = new double [pd_ptr->n_pair];
   fcst_rank = new double [pd_ptr->n_pair];
   obs_rank  = new double [pd_ptr->n_pair];
   
   // Initialize the sums
   count   = 0;
   f_sum   = o_sum       = ff_sum     = oo_sum = fo_sum = 0.0;
   err_sum = abs_err_sum = err_sq_sum = 0.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];
      
      // Skip bad data values
      if(f == bad_data_double || o == bad_data_double) continue;
      
      // 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;
   } // end for i
   
   if(count == 0) {
      cerr << "\n\nERROR: compute_cnt() -> "
           << "count is zero!\n\n" << flush;
      exit(1);
   }
   
   // Store the sample size
   cnt_info.n = count;
   
   // Compute the forecast mean
   cnt_info.fbar.v = f_sum/count;

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

   // Compute the observation mean
   cnt_info.obar.v = o_sum/count;

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

   // Compute the f*o mean
   cnt_info.fobar = fo_sum/count;
   
   // Compute the forecast squared mean
   cnt_info.ffbar = ff_sum/count;
   
   // Compute the observation squared mean
   cnt_info.oobar = oo_sum/count;
   
   // Comptue the bias
   cnt_info.bias = cnt_info.fbar.v - cnt_info.obar.v;
   
   // Compute multiplicative bias
   if(cnt_info.obar.v == 0) 
      cnt_info.mbias = bad_data_double;
   else
      cnt_info.mbias = cnt_info.fbar.v/cnt_info.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.pr_corr.v = bad_data_double; 
   else           cnt_info.pr_corr.v = ((count*fo_sum) - (f_sum*o_sum))/den;
   
   // Compute the Kendall Tau and Spearman Rank correlation coefficients
      
   // Compute the rank of the raw fcst and obs data values
   cnt_info.n_ranks    = pd_ptr->n_pair;
   cnt_info.frank_ties = rank(pd_ptr->fdata, fcst_rank, pd_ptr->n_pair);
   cnt_info.orank_ties = rank(pd_ptr->odata, obs_rank, pd_ptr->n_pair);
   
   // Compute the sums for the ranks for use in computing Spearman's
   // Rank correlation coefficient
   f_sum = o_sum = ff_sum = oo_sum = fo_sum = 0.0;
   for(i=0; i<pd_ptr->n_pair; i++) {
   
      f = fcst_rank[i];
      o = obs_rank[i];

      f_sum  += f;
      o_sum  += o;
      ff_sum += f*f;
      oo_sum += o*o;
      fo_sum += f*o;
   } // end for i
   
   // 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.sp_corr = bad_data_double; 
   else           cnt_info.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<pd_ptr->n_pair; j++) {
      for(k=j+1; k<pd_ptr->n_pair; 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.kt_corr = bad_data_double; 
   else           cnt_info.kt_corr = (concordant - discordant)/den;

   // Compute the percntiles of the error
   sort(err_array, count);
   cnt_info.e10 = percentile(err_array, count, 0.10);
   cnt_info.e25 = percentile(err_array, count, 0.25);
   cnt_info.e50 = percentile(err_array, count, 0.50);
   cnt_info.e75 = percentile(err_array, count, 0.75);
   cnt_info.e90 = percentile(err_array, count, 0.90);

   // Compute mean error and standard deviation of the mean error
   cnt_info.me.v = err_sum/count;
   cnt_info.estdev.v = compute_stdev(err_sum, err_sq_sum, count);
      
   // Compute mean absolute error
   cnt_info.mae = abs_err_sum/count;
      
   // Compute the mean squared error
   cnt_info.mse = err_sq_sum/count;
   
   // Compute the bias corrected mean squared error (decomposition of MSE)
   f = cnt_info.fbar.v;
   o = cnt_info.obar.v;
   cnt_info.bcmse = cnt_info.mse - (f-o)*(f-o); 
   
   // Compute root mean squared error
   cnt_info.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(i=0; i<n_ci_alpha; i++) {

      // Store the alpha value in the CNTInfo object
      cnt_info.alpha[i] = conf.ci_alpha(i).dval();
      
      // Check for the degenerate case
      if(count <= 1) {
         cnt_info.fbar.v_cl[i]    = cnt_info.fbar.v_cu[i]    = bad_data_double;
         cnt_info.fstdev.v_cl[i]  = cnt_info.fstdev.v_cu[i]  = bad_data_double;
         cnt_info.obar.v_cl[i]    = cnt_info.obar.v_cu[i]    = bad_data_double;
         cnt_info.ostdev.v_cl[i]  = cnt_info.ostdev.v_cu[i]  = bad_data_double;
         cnt_info.pr_corr.v_cl[i] = cnt_info.pr_corr.v_cu[i] = bad_data_double;
         cnt_info.me.v_cl[i]      = cnt_info.me.v_cu[i]      = bad_data_double;
         cnt_info.estdev.v_cl[i]  = cnt_info.estdev.v_cu[i]  = 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.alpha[i]/2.0, 0.0, 1.0);
         cv_normal_u = normal_cdf_inv(1.0 - (cnt_info.alpha[i]/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.alpha[i]/2.0, count - 1);
         cv_normal_u = students_t_cdf_inv(1.0 - (cnt_info.alpha[i]/2.0), count - 1);
      }

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

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

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

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

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

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

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

      // Compute confidence interval for the error standard deviation
      cnt_info.estdev.v_cl[i] = sqrt((count - 1)*cnt_info.estdev.v*cnt_info.estdev.v/cv_chi2_u);
      cnt_info.estdev.v_cu[i] = sqrt((count - 1)*cnt_info.estdev.v*cnt_info.estdev.v/cv_chi2_l);
   } // end for i  
   
   // Deallocate memory
   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 compute_sl1l2(SL1L2Info &s_info, int i_gc, PairData *pd_ptr) {
   int i, status;
   double f, c, o;
   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) {
   
      strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, conf.ncep_defaults().ival(), status));
      
      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;
   }
      
   // 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];
      c = pd_ptr->cdata[i];
      o = pd_ptr->odata[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, status;
   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) {
   
      strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_ugrd].gc.code, conf.ncep_defaults().ival(), status));
      strcpy(abbr_str_2, get_grib_code_abbr(gc_pd[i_vgrd].gc.code, conf.ncep_defaults().ival(), status));
      
      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);
   }
   
   // 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_cts(CTSInfo &cts_info, int i_gc, PairData *pd_ptr) {
   int j, mon, day, yr, hr, min, sec, status;
   int l_hr, l_min, l_sec;
   char line[max_line_len], fmt_str[max_str_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];

   // Compute a string for the lead time and valid time
   unix_to_mdyhms(gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(), 
                  mon, day, yr, hr, min, sec);
   sec_to_hms(gc_pd[i_gc].fcst_wd_ptr[0]->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);

   // Type of thresholding performed
   if(cts_info.cts_thresh_ind >= 0
   && cts_info.cts_thresh_ind < n_thresh_type) {
      sprintf(thresh_str, "%s%.3f", 
         thresh_type_str[cts_info.cts_thresh_ind], 
         cts_info.cts_thresh);
   }
   else {
      cerr << "\n\nERROR: write_cts() -> "
           << "unexpected threshold indicator value of " 
           << cts_info.cts_thresh_ind << "\".\n\n" << flush;
      exit(1);
   }
   
   // Grib Field abbreviation
   strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, 
          conf.ncep_defaults().ival(), status));

   // Interpolation method and number of points
   strcpy(interp_mthd_str, mthd_str[pd_ptr->interp_mthd]);
   sprintf(interp_pnts_str, "%i", 
           pd_ptr->interp_wdth*pd_ptr->interp_wdth);

   // 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.cts.n());
      dbl_to_str(cts_info.cts.f_rate(status), var_str[1]);
      dbl_to_str(cts_info.cts.h_rate(status), var_str[2]);
      dbl_to_str(cts_info.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
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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.cts.n());
      sprintf(var_str[1], "%i", cts_info.cts.fy_oy());
      sprintf(var_str[2], "%i", cts_info.cts.fy_on());
      sprintf(var_str[3], "%i", cts_info.cts.fn_oy());
      sprintf(var_str[4], "%i", cts_info.cts.fn_on());
      
      // Setup the format string
      ctc_fmt_str(fmt_str);

      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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.cts.n());
      dbl_to_str(cts_info.cts.fy_oy_tp(status), var_str[1]);
      dbl_to_str(cts_info.cts.fy_on_tp(status), var_str[2]);
      dbl_to_str(cts_info.cts.fn_oy_tp(status), var_str[3]);
      dbl_to_str(cts_info.cts.fn_on_tp(status), var_str[4]);
      dbl_to_str(cts_info.cts.fy_tp(status), var_str[5]);
      dbl_to_str(cts_info.cts.fn_tp(status), var_str[6]);
      dbl_to_str(cts_info.cts.oy_tp(status), var_str[7]);
      dbl_to_str(cts_info.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
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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.cts.n());
      dbl_to_str(cts_info.cts.fy_oy_fp(status), var_str[1]);
      dbl_to_str(cts_info.cts.fy_on_fp(status), var_str[2]);
      dbl_to_str(cts_info.cts.fn_oy_fp(status), var_str[3]);
      dbl_to_str(cts_info.cts.fn_on_fp(status), var_str[4]);
      sprintf(var_str[5], "%i", cts_info.cts.fy());
      sprintf(var_str[6], "%i", cts_info.cts.fn());

      // Setup the format string
      cfp_fmt_str(fmt_str);

      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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.cts.n());
      dbl_to_str(cts_info.cts.fy_oy_op(status), var_str[1]);
      dbl_to_str(cts_info.cts.fy_on_op(status), var_str[2]);
      dbl_to_str(cts_info.cts.fn_oy_op(status), var_str[3]);
      dbl_to_str(cts_info.cts.fn_on_op(status), var_str[4]);
      sprintf(var_str[5], "%i", cts_info.cts.oy());
      sprintf(var_str[6], "%i", cts_info.cts.on());

      // Setup the format string
      cop_fmt_str(fmt_str);

      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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.n_alpha; j++) {

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

         // Convert int and double values to strings
         sprintf(var_str[0], "%i", cts_info.cts.n());
         dbl_to_str(cts_info.baser.v, var_str[1]);
         dbl_to_str(cts_info.baser.v_cl[j], var_str[2]);
         dbl_to_str(cts_info.baser.v_cu[j], var_str[3]);
         dbl_to_str(cts_info.fmean.v, var_str[4]);
         dbl_to_str(cts_info.fmean.v_cl[j], var_str[5]);
         dbl_to_str(cts_info.fmean.v_cu[j], var_str[6]);
         dbl_to_str(cts_info.acc.v, var_str[7]);
         dbl_to_str(cts_info.acc.v_cl[j], var_str[8]);
         dbl_to_str(cts_info.acc.v_cu[j], var_str[9]);
         dbl_to_str(cts_info.fbias.v, var_str[10]);
         dbl_to_str(cts_info.pody.v, var_str[11]);
         dbl_to_str(cts_info.pody.v_cl[j], var_str[12]);
         dbl_to_str(cts_info.pody.v_cu[j], var_str[13]);
         dbl_to_str(cts_info.podn.v, var_str[14]);
         dbl_to_str(cts_info.podn.v_cl[j], var_str[15]);
         dbl_to_str(cts_info.podn.v_cu[j], var_str[16]);
         dbl_to_str(cts_info.pofd.v, var_str[17]);
         dbl_to_str(cts_info.pofd.v_cl[j], var_str[18]);
         dbl_to_str(cts_info.pofd.v_cu[j], var_str[19]);
         dbl_to_str(cts_info.far.v, var_str[20]);
         dbl_to_str(cts_info.far.v_cl[j], var_str[21]);
         dbl_to_str(cts_info.far.v_cu[j], var_str[22]);
         dbl_to_str(cts_info.csi.v, var_str[23]);
         dbl_to_str(cts_info.csi.v_cl[j], var_str[24]);
         dbl_to_str(cts_info.csi.v_cu[j], var_str[25]);
         dbl_to_str(cts_info.gss.v, var_str[26]);
         dbl_to_str(cts_info.hk.v, var_str[27]);
         dbl_to_str(cts_info.hk.v_cl[j], var_str[28]);
         dbl_to_str(cts_info.hk.v_cu[j], var_str[29]);
         dbl_to_str(cts_info.hss.v, var_str[30]);
         dbl_to_str(cts_info.odds.v, var_str[31]);
         dbl_to_str(cts_info.odds.v_cl[j], var_str[32]);
         dbl_to_str(cts_info.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
            conf.model().sval(),                   // Model Name
            lead_str,                              // Lead Time HHMMSS
            valid_str,                             // Valid Time YYYYMMDD_HHMMSS
            pd_ptr->msg_typ,                       // Verifying Observation Type
            pd_ptr->mask_name,                     // Verification masking region
            type_str,                              // Threshold applied
            abbr_str,                              // Grib field name
            gc_pd[i_gc].gc.lvl_str,                // 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

   return;
}

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

void write_cnt(CNTInfo &cnt_info, int i_gc, PairData *pd_ptr) {
   int j, mon, day, yr, hr, min, sec, status;
   int l_hr, l_min, l_sec;
   char line[max_line_len], fmt_str[max_str_len];
   char lead_str[max_str_len], valid_str[max_str_len], type_str[max_str_len];
   char var_str[max_num_vars][var_str_len];
   char interp_mthd_str[max_str_len], interp_pnts_str[max_str_len];

   // Compute a string for the lead time and valid time
   unix_to_mdyhms(gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(), 
                  mon, day, yr, hr, min, sec);
   sec_to_hms(gc_pd[i_gc].fcst_wd_ptr[0]->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);   
   
   // Grib Field abbreviation
   strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, 
          conf.ncep_defaults().ival(), status));

   // Interpolation method and number of points
   strcpy(interp_mthd_str, mthd_str[pd_ptr->interp_mthd]);
   sprintf(interp_pnts_str, "%i", 
           pd_ptr->interp_wdth*pd_ptr->interp_wdth);

   // Dump out a VSDB record for each Continuous Statistics object
   for(j=0; j<n_ci_alpha; j++) {
   
      // 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) {
      
         sprintf(type_str, "CNT/%.3f", cnt_info.alpha[j]);

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

         // Setup the format string
         cnt_fmt_str(fmt_str);

         // Write the line
         sprintf(line, fmt_str,
            met_version,                           // MET Version
            conf.model().sval(),                   // Model Name
            lead_str,                              // Lead Time HHMMSS
            valid_str,                             // Valid Time YYYYMMDD_HHMMSS
            pd_ptr->msg_typ,                       // Verifying Observation Type
            pd_ptr->mask_name,                     // Verification masking region
            type_str,                              // Threshold applied
            abbr_str,                              // Grib field name
            gc_pd[i_gc].gc.lvl_str,                // 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 if
   } // end for i

   return;
}

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

void write_sl1l2(SL1L2Info &s_info, int i_gc, PairData *pd_ptr) {
   int mon, day, yr, hr, min, sec, status;
   int l_hr, l_min, l_sec;
   char line[max_line_len], fmt_str[max_str_len];
   char lead_str[max_str_len], valid_str[max_str_len], type_str[max_str_len];
   char var_str[max_num_vars][var_str_len];
   char interp_mthd_str[max_str_len], interp_pnts_str[max_str_len];

   // Compute a string for the lead time and valid time
   unix_to_mdyhms(gc_pd[i_gc].fcst_wd_ptr[0]->get_valid_time(), 
                  mon, day, yr, hr, min, sec);
   sec_to_hms(gc_pd[i_gc].fcst_wd_ptr[0]->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);   
   
   // Grib Field abbreviation
   strcpy(abbr_str, get_grib_code_abbr(gc_pd[i_gc].gc.code, 
          conf.ncep_defaults().ival(), status));

   // Interpolation method and number of points
   strcpy(interp_mthd_str, mthd_str[pd_ptr->interp_mthd]);
   sprintf(interp_pnts_str, "%i", 
           pd_ptr->interp_wdth*pd_ptr->interp_wdth);   
      
   // 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 && s_info.scount > 0) {
   
      strcpy(type_str, "SL1L2");

      // Convert double values to strings
      sprintf(var_str[0], "%i", s_info.scount);
      dbl_to_str(s_info.fbar,   var_str[1]);
      dbl_to_str(s_info.obar,   var_str[2]);
      dbl_to_str(s_info.fobar,  var_str[3]);
      dbl_to_str(s_info.ffbar,  var_str[4]);
      dbl_to_str(s_info.oobar,  var_str[5]);

      // Setup the format string
      sl1l2_fmt_str(fmt_str);
      
      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // 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";
   }
   
   // Scalar Anomaly L1L2 Line Type (SAL1L2)
   // Dump out the SL1L2 VSDB line:
   //    TOTAL, FABAR, OABAR, FOABAR, FFABAR, OOABAR,
   //    INTERP_MTHD, INTERP_PNTS
   if(conf.output_flag(i_sal1l2).ival() != 0 && s_info.sacount > 0) {
   
      strcpy(type_str, "SAL1L2");

      // Convert double values to strings
      sprintf(var_str[0], "%i",  s_info.sacount);
      dbl_to_str(s_info.fabar,   var_str[1]);
      dbl_to_str(s_info.oabar,   var_str[2]);
      dbl_to_str(s_info.foabar,  var_str[3]);
      dbl_to_str(s_info.ffabar,  var_str[4]);
      dbl_to_str(s_info.ooabar,  var_str[5]);

      // Setup the format string
      sal1l2_fmt_str(fmt_str);

      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         pd_ptr->msg_typ,                   // Verifying Observation Type
         pd_ptr->mask_name,                 // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_gc].gc.lvl_str,            // Grib level name
         var_str[0],                        // Total anomaly count
         var_str[1],                        // fabar
         var_str[2],                        // oabar
         var_str[3],                        // foabar
         var_str[4],                        // ffabar
         var_str[5],                        // ooabar
         interp_mthd_str,                   // Interpolation method 
         interp_pnts_str                    // Interpolation points
      );
      *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) {
   int mon, day, yr, hr, min, sec, status;
   int l_hr, l_min, l_sec;
   char line[max_line_len], fmt_str[max_str_len];
   char lead_str[max_str_len], valid_str[max_str_len], type_str[max_str_len];
   char var_str[max_num_vars][var_str_len];
   char interp_mthd_str[max_str_len], interp_pnts_str[max_str_len];
   char tmp_str[max_str_len], tmp_str_2[max_str_len];

   // Compute a string for the lead time and valid time
   unix_to_mdyhms(gc_pd[i_ugrd].fcst_wd_ptr[0]->get_valid_time(), 
                  mon, day, yr, hr, min, sec);
   sec_to_hms(gc_pd[i_ugrd].fcst_wd_ptr[0]->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);   
   
   // Grib Field abbreviation
   strcpy(tmp_str, 
          get_grib_code_abbr(gc_pd[i_ugrd].gc.code, conf.ncep_defaults().ival(), status));
   strcpy(tmp_str_2, 
          get_grib_code_abbr(gc_pd[i_vgrd].gc.code, conf.ncep_defaults().ival(), status));
   sprintf(abbr_str, "%s_%s", tmp_str, tmp_str_2);

   // Interpolation method and number of points
   strcpy(interp_mthd_str, mthd_str[ugrd_pd_ptr->interp_mthd]);
   sprintf(interp_pnts_str, "%i", 
           ugrd_pd_ptr->interp_wdth*ugrd_pd_ptr->interp_wdth);   
      
   // Vector L1L2 Line Type (VL1L2)
   // Dump out the VL1L2 VSDB line:
   //    TOTAL, UFBAR, VFBAR, UOBAR, VOBAR, UVFOBAR, UVFFBAR, UVOOBAR,
   //    INTERP_MTHD, INTERP_PNTS
   if(conf.output_flag(i_vl1l2).ival() != 0 && v_info.vcount > 0) {
      
      strcpy(type_str, "VL1L2");

      // Convert double values to strings
      sprintf(var_str[0], "%i",  v_info.vcount);
      dbl_to_str(v_info.ufbar,   var_str[1]);
      dbl_to_str(v_info.vfbar,   var_str[2]);
      dbl_to_str(v_info.uobar,   var_str[3]);
      dbl_to_str(v_info.vobar,   var_str[4]);
      dbl_to_str(v_info.uvfobar, var_str[5]);
      dbl_to_str(v_info.uvffbar, var_str[6]);
      dbl_to_str(v_info.uvoobar, var_str[7]);
      
      // Setup the format string
      vl1l2_fmt_str(fmt_str);
      
      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         ugrd_pd_ptr->msg_typ,              // Verifying Observation Type
         ugrd_pd_ptr->mask_name,            // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_ugrd].gc.lvl_str,          // Grib level name
         var_str[0],                        // Total count
         var_str[1],                        // ufbar
         var_str[2],                        // vfbar
         var_str[3],                        // uobar
         var_str[4],                        // vobar
         var_str[5],                        // uvfobar
         var_str[6],                        // uvffbar
         var_str[7],                        // uvoobar
         interp_mthd_str,                   // Interpolation method 
         interp_pnts_str                    // Interpolation points
      );
      *vsdb_out << line << "\n";

      if(conf.output_flag(i_vl1l2).ival() >= 2) *vl1l2_out << line << "\n";
   }
   
   // Vector Anomaly L1L2 Line Type (VAL1L2)
   // Dump out the VAL1L2 VSDB line:
   //    TOTAL, UFABAR, VFABAR, UOABAR, VOABAR, UVFOABAR, UVFFABAR, UVOOABAR,
   //    INTERP_MTHD, INTERP_PNTS
   if(conf.output_flag(i_val1l2).ival() != 0 && v_info.vacount > 0) {
   
      strcpy(type_str, "VAL1L2");

      // Convert double values to strings
      sprintf(var_str[0], "%i",   v_info.vacount);
      dbl_to_str(v_info.ufabar,   var_str[1]);
      dbl_to_str(v_info.vfabar,   var_str[2]);
      dbl_to_str(v_info.uoabar,   var_str[3]);
      dbl_to_str(v_info.voabar,   var_str[4]);
      dbl_to_str(v_info.uvfoabar, var_str[5]);
      dbl_to_str(v_info.uvffabar, var_str[6]);
      dbl_to_str(v_info.uvooabar, var_str[7]);
      
      // Setup the format string
      val1l2_fmt_str(fmt_str);
      
      // Write the line
      sprintf(line, fmt_str,
         met_version,                       // MET Version
         conf.model().sval(),               // Model Name
         lead_str,                          // Lead Time HHMMSS
         valid_str,                         // Valid Time YYYYMMDD_HHMMSS
         ugrd_pd_ptr->msg_typ,              // Verifying Observation Type
         ugrd_pd_ptr->mask_name,            // Verification masking region
         type_str,                          // Threshold applied
         abbr_str,                          // Grib field name
         gc_pd[i_ugrd].gc.lvl_str,          // Grib level name
         var_str[0],                        // Total anomaly count
         var_str[1],                        // ufabar
         var_str[2],                        // vfabar
         var_str[3],                        // uoabar
         var_str[4],                        // voabar
         var_str[5],                        // uvfoabar
         var_str[6],                        // uvffabar
         var_str[7],                        // uvooabar
         interp_mthd_str,                   // Interpolation method 
         interp_pnts_str                    // Interpolation points
      );  
      *vsdb_out << line << "\n";

      if(conf.output_flag(i_val1l2).ival() >= 2) *val1l2_out << line << "\n";
   }
   
   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;
   }

   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[-outdir path]\n"
        << "\t[-v level]\n\n"

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

        << "\t\t\"obs_file\" is an observation file in netCDF format "
        << "(output of pb2nc) 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\" provides a Grib file containing "
        << "climatological values on the same grid as the forecast file\n"
        << "\t\t       to be used when computing scalar and vector "
        << "anomaly measures using the contents of \"climo_file\".\n"
        << "\t\t       If not provided, scalar and vector anomaly values "
        << "will not be computed (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;
}

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