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

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

using namespace std;

#include <iostream>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>

#include "vx_met_util/met_stats.h"
#include "vx_met_util/compute_ci.h"
#include "vx_met_util/constants.h"
#include "vx_util/vx_util.h"
#include "vx_grib_classes/grib_strings.h"
#include "vx_wrfdata/vx_wrfdata.h"

////////////////////////////////////////////////////////////////////////
//
// Code for class ThreshInfo
//
////////////////////////////////////////////////////////////////////////

ThreshInfo::ThreshInfo() {
   init_from_scratch();
}

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

ThreshInfo::~ThreshInfo() {
   clear();
}

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

ThreshInfo::ThreshInfo(const ThreshInfo &c) {
   
   init_from_scratch();

   assign(c);
}

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

ThreshInfo & ThreshInfo::operator=(const ThreshInfo &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

void ThreshInfo::init_from_scratch() {

   thresh     = (double *) 0;
   thresh_ind = (int *) 0;   

   clear();

   return;
}

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

void ThreshInfo::clear() {

   if(thresh)     { delete [] thresh;     thresh = (double *) 0; }
   if(thresh_ind) { delete [] thresh_ind; thresh_ind = (int *) 0; }
   
   n_thresh = 0;
   
   return;
}

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

void ThreshInfo::assign(const ThreshInfo &ti) {
   int i;

   clear();
   
   set_n_thresh(ti.n_thresh);
   
   for(i=0; i<ti.n_thresh; i++) {
      set_thresh(i, ti.thresh[i]);
      set_thresh_ind(i, ti.thresh_ind[i]);
   }

   return;
}

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

void ThreshInfo::set_n_thresh(int n) {

   n_thresh = n;

   thresh     = new double [n_thresh];
   thresh_ind = new int [n_thresh];
   
   if(!thresh || !thresh_ind) {
      cerr << "\n\nERROR: ThreshInfo::set_n_thresh() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }

   return;
}

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

void ThreshInfo::set_thresh(int i, double t) {

   thresh[i] = t;

   return;
}

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

void ThreshInfo::set_thresh_ind(int i, int ind) {

   thresh_ind[i] = ind;

   return;
}

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

void ThreshInfo::parse_thresh_str(const char *thresh_str) {
   char tmp_str[max_str_len], a, b;
   char *c = (char *) 0;
   int i, n;

   // Compute the number of tokens in the string based on " "
   n = num_tokens(thresh_str, " ");
   
   // Check for no tokens in string
   if(n == 0) return;
   
   // Allocate space for all of the threshold values
   set_n_thresh(n);   
   
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, thresh_str);

   // Tokenize the string and store the threshold values
   c = strtok(tmp_str, " ");
   
   // Check that the threshold string is more than 2 characters long
   if(strlen(c) < 3) {
      cerr << "\n\nERROR: ThreshInfo::parse_thresh_str() -> "
           << "each threshold value must be preceeded by one of "
           << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
           << "thresholding to be performed \"" << thresh_str
           << "\".\n\n" << flush;
      exit(1);
   }
   
   a = tolower(c[0]);
   b = tolower(c[1]);
   
   if     (a == 'l' && b == 't') thresh_ind[0] = thresh_lt; // less than
   else if(a == 'l' && b == 'e') thresh_ind[0] = thresh_le; // less than or equal to
   else if(a == 'e' && b == 'q') thresh_ind[0] = thresh_eq; // equal to
   else if(a == 'n' && b == 'e') thresh_ind[0] = thresh_ne; // not equal to
   else if(a == 'g' && b == 't') thresh_ind[0] = thresh_gt; // greater than
   else if(a == 'g' && b == 'e') thresh_ind[0] = thresh_ge; // greater than or equal to
   else {
      cerr << "\n\nERROR: ThreshInfo::parse_thresh_str() -> "
           << "each threshold value must be preceeded by one of "
           << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
           << "thresholding to be performed \"" << thresh_str
           << "\".\n\n" << flush;
      exit(1);   
   }
   c += 2;
   thresh[0] = atof(c);
   
   // Parse remaining tokens
   for(i=1; i<n; i++) {
      
      c = strtok(0, " ");
      
      // Check that the threshold string is more than 2 characters long
      if(strlen(c) < 3) {
         cerr << "\n\nERROR: ThreshInfo::parse_thresh_str() -> "
              << "each threshold value must be preceeded by one of "
              << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
              << "thresholding to be performed \"" << thresh_str
              << "\".\n\n" << flush;
         exit(1);
      }
      
      a = tolower(c[0]);
      b = tolower(c[1]);
   
      if     (a == 'l' && b == 't') thresh_ind[i] = thresh_lt; // less than
      else if(a == 'l' && b == 'e') thresh_ind[i] = thresh_le; // less than or equal to
      else if(a == 'e' && b == 'q') thresh_ind[i] = thresh_eq; // equal to
      else if(a == 'n' && b == 'e') thresh_ind[i] = thresh_ne; // not equal to
      else if(a == 'g' && b == 't') thresh_ind[i] = thresh_gt; // greater than
      else if(a == 'g' && b == 'e') thresh_ind[i] = thresh_ge; // greater than or equal to
      else {
         cerr << "\n\nERROR: ThreshInfo::parse_thresh_str() -> "
              << "each threshold value must be preceeded by one of "
              << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
              << "thresholding to be performed \"" << thresh_str
              << "\".\n\n" << flush;
         exit(1);   
      }
      c += 2;
      thresh[i] = atof(c);
   }
   
   return;  
}

////////////////////////////////////////////////////////////////////////
//
// Code for class GCInfo
//
////////////////////////////////////////////////////////////////////////

GCInfo::GCInfo() {
   init_from_scratch();
}

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

GCInfo::~GCInfo() {
   clear();
}

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

GCInfo::GCInfo(const GCInfo &c) {
   
   init_from_scratch();

   assign(c);
}

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

GCInfo & GCInfo::operator=(const GCInfo &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

int GCInfo::operator==(const GCInfo &c) {
   int match = 0;

   if(code        == c.code &&
      lvl_type    == c.lvl_type &&
      lvl_1       == c.lvl_1 &&
      lvl_2       == c.lvl_2 &&
      vector_flag == c.vector_flag) match = 1;

   return(match);
}


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

void GCInfo::init_from_scratch() {

   lvl_str = (char *) 0;

   clear();

   return;
}

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

void GCInfo::clear() {

   code        = 0;
   lvl_type    = NoLevel;
   lvl_1       = 0;
   lvl_2       = 0;
   vector_flag = 0;
   
   if(lvl_str) { delete [] lvl_str; lvl_str = (char *) 0; }

   return;
}

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

void GCInfo::assign(const GCInfo &c) {

   clear();
   
   code        = c.code;
   lvl_type    = c.lvl_type;
   lvl_1       = c.lvl_1;
   lvl_2       = c.lvl_2;
   vector_flag = c.vector_flag;
   
   set_lvl_str(c.lvl_str);

   return;
}

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

void GCInfo::set_lvl_str(const char *c) {

   lvl_str = new char [strlen(c)+1];
   
   if(!lvl_str) {
      cerr << "\n\nERROR: GCInfo::set_lvl_str() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }
   
   strcpy(lvl_str, c);

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class PairData
//
////////////////////////////////////////////////////////////////////////

PairData::PairData() {
   init_from_scratch();
}

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

PairData::~PairData() {
   clear();
}

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

PairData::PairData(const PairData &pd) {
   
   init_from_scratch();

   assign(pd);
}

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

PairData & PairData::operator=(const PairData &pd) {

   if(this == &pd) return(*this);

   assign(pd);

   return(*this);
}

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

void PairData::init_from_scratch() {

   mask_name    = (char *) 0;
   mask_wd_ptr  = (WrfData *) 0;
   
   msg_typ      = (char *) 0;
   
   interp_mthd  = -1;
   interp_wdth  = -1;
   
   fdata        = (double *) 0;
   cdata        = (double *) 0;
   odata        = (double *) 0;

   n_pair       = 0;
   n_alloc      = 0;
   
   clear();
   
   return;
}

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

void PairData::clear() {

   if(mask_name) { delete [] mask_name; mask_name = (char *) 0; }
   
   mask_wd_ptr = (WrfData *) 0;  // Not allocated
   
   if(msg_typ)   { delete [] msg_typ;   msg_typ = (char *) 0; }
   
   interp_mthd = -1;
   interp_wdth = -1;
   
   if(fdata)     { delete [] fdata;     fdata = (double *) 0; }
   if(cdata)     { delete [] cdata;     cdata = (double *) 0; }
   if(odata)     { delete [] odata;     odata = (double *) 0; }
   
   n_pair  = 0;
   n_alloc = 0;
   
   return;
}

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

void PairData::assign(const PairData &pd) {
   int i;

   clear();
   
   set_mask_name(pd.mask_name);
   set_mask_wd_ptr(pd.mask_wd_ptr);
   
   set_msg_typ(pd.msg_typ);
   
   set_interp_mthd(pd.interp_mthd);
   set_interp_wdth(pd.interp_wdth);
   
   for(i=0; i<pd.n_pair; i++) {
      add_pair(pd.fdata[i], pd.cdata[i], pd.odata[i]);
   }

   return;
}

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

void PairData::set_mask_name(const char *c) {
   
   mask_name = new char [strlen(c)+1];
   
   if(!mask_name) {
      cerr << "\n\nERROR: PairData::set_mask_name() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }
   
   strcpy(mask_name, c);

   return;
}

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

void PairData::set_mask_wd_ptr(WrfData *wd_ptr) {

   mask_wd_ptr = wd_ptr;
   
   return;
}

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

void PairData::set_msg_typ(const char *c) {
   
   msg_typ = new char [strlen(c)+1];
   
   if(!msg_typ) {
      cerr << "\n\nERROR: PairData::set_msg_typ() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }
   
   strcpy(msg_typ, c);

   return;
}

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

void PairData::set_interp_mthd(const char *str) {
   int i, n;

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

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

void PairData::set_interp_mthd(int m) {

   interp_mthd = m;
   
   return;
}

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

void PairData::set_interp_wdth(int n) {

   interp_wdth = n;
   
   return;
}

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

void PairData::add_pair(double fval, double cval, double oval) {
   int i;
   double *f_tmp, *c_tmp, *o_tmp;
   
   // If the number of pairs exceeds the amount of allocated space
   // reallocate more space
   if(n_pair >= n_alloc) {
      
      f_tmp = new double [n_alloc + pair_data_alloc_jump];
      c_tmp = new double [n_alloc + pair_data_alloc_jump];
      o_tmp = new double [n_alloc + pair_data_alloc_jump];
      
      if(!f_tmp || !c_tmp || !o_tmp) {
         cerr << "\n\nERROR: PairData::add_pair() -> "
              << "Memory allocation error!\n\n" << flush;
         exit(1);
      }
      
      // Copy over the existing data
      for(i=0; i<n_pair; i++) {
         f_tmp[i] = fdata[i];
         c_tmp[i] = cdata[i];
         o_tmp[i] = odata[i];
      }
      
      // Delete the old fdata, cdata, and odata arrays
      if(fdata) { delete [] fdata; fdata = (double *) 0; }
      if(cdata) { delete [] cdata; cdata = (double *) 0; }
      if(odata) { delete [] odata; odata = (double *) 0; }
      
      // Point the old fdata, cdata, and odata arrays to the newly 
      // allocated space
      fdata = f_tmp;
      cdata = c_tmp;
      odata = o_tmp;
   
      // Increment the allocation count
      n_alloc += pair_data_alloc_jump;
   }
   
   // Save the pair data
   fdata[n_pair] = fval;
   cdata[n_pair] = cval;
   odata[n_pair] = oval;
   
   // Increment the number of pairs
   n_pair += 1;
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class GCPairData
//
////////////////////////////////////////////////////////////////////////

GCPairData::GCPairData() {
   init_from_scratch();
}

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

GCPairData::~GCPairData() {
   clear();
}

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

GCPairData::GCPairData(const GCPairData &gc_pd) {
   
   init_from_scratch();

   assign(gc_pd);
}

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

GCPairData & GCPairData::operator=(const GCPairData &gc_pd) {

   if(this == &gc_pd) return(*this);

   assign(gc_pd);

   return(*this);
}

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

void GCPairData::init_from_scratch() {
   
   interp_threshold = 0;

   n_fcst           = 0;   
   fcst_prs_lvl     = (int *) 0;
   fcst_wd_ptr      = (WrfData **) 0;
   
   n_climo          = 0;
   climo_prs_lvl    = (int *) 0;
   climo_wd_ptr     = (WrfData **) 0;
   
   n_msg_typ        = 0;
   n_mask           = 0;
   n_interp         = 0;
   
   pd               = (PairData ***) 0;
   
   clear();
   
   return;
}

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

void GCPairData::clear() {
   int i, j, k;

   gc.clear();
   
   interp_threshold = 0;
   
   if(fcst_prs_lvl)  { delete [] fcst_prs_lvl;  fcst_prs_lvl  = (int *) 0; }
   if(climo_prs_lvl) { delete [] climo_prs_lvl; climo_prs_lvl = (int *) 0; }
   
   for(i=0; i<n_fcst; i++)  fcst_wd_ptr[i]  = (WrfData *) 0;
   for(i=0; i<n_climo; i++) climo_wd_ptr[i] = (WrfData *) 0;
   
   fcst_wd_ptr  = (WrfData **) 0;
   climo_wd_ptr = (WrfData **) 0;

   for(i=0; i<n_msg_typ; i++) {
      for(j=0; j<n_mask; j++) {
         for(k=0; k<n_interp; k++) {     
            pd[i][j][k].clear();
         }
      }
   }
   
   n_fcst    = 0;
   n_climo   = 0;
   n_msg_typ = 0;
   n_mask    = 0;
   n_interp  = 0;
   
   return;
}

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

void GCPairData::assign(const GCPairData &gc_pd) {
   int i, j, k;

   clear();
   
   set_gc(gc_pd.gc);
   
   interp_threshold = gc_pd.interp_threshold;
   
   set_n_fcst(gc_pd.n_fcst);
   for(i=0; i<gc_pd.n_fcst; i++) {
      set_fcst_prs_lvl(i, fcst_prs_lvl[i]);
      set_fcst_wd_ptr(i, fcst_wd_ptr[i]);
   }
   
   set_n_climo(gc_pd.n_climo);
   for(i=0; i<gc_pd.n_climo; i++) {
      set_climo_prs_lvl(i, climo_prs_lvl[i]);
      set_climo_wd_ptr(i, climo_wd_ptr[i]);
   }
   
   set_pd_size(gc_pd.n_msg_typ, gc_pd.n_mask, gc_pd.n_interp);
      
   for(i=0; i<gc_pd.n_msg_typ; i++) {
      for(j=0; j<gc_pd.n_mask; j++) {
         for(k=0; k<gc_pd.n_interp; k++) {

            pd[i][j][k] = gc_pd.pd[i][j][k];
         }
      }
   }
   
   return;
}
     
////////////////////////////////////////////////////////////////////////

void GCPairData::set_gc(const GCInfo &gc_new) {

   gc = gc_new;
   
   return;
}

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

void GCPairData::set_interp_threshold(double t) {

   interp_threshold = t;
   
   return;
}

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

void GCPairData::set_n_fcst(int n) {
   int i;

   n_fcst = n;
   
   fcst_prs_lvl = new int       [n_fcst];
   fcst_wd_ptr  = new WrfData * [n_fcst];
   
   if(!fcst_prs_lvl || !fcst_wd_ptr ) {
      cerr << "\n\nERROR: GCPairData::set_n_fcst() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }
   
   for(i=0; i<n_fcst; i++) {
      fcst_prs_lvl[i] = bad_data_int;
      fcst_wd_ptr[i]  = (WrfData *) 0;
   }
   
   return;
}

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

void GCPairData::set_fcst_prs_lvl(int i, int lvl) {

   fcst_prs_lvl[i] = lvl;
   
   return;
}

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

void GCPairData::set_fcst_wd_ptr(int i, WrfData *wd_ptr) {

   fcst_wd_ptr[i] = wd_ptr;
   
   return;
}

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

void GCPairData::set_n_climo(int n) {
   int i;

   n_climo = n;
   
   climo_prs_lvl = new int       [n_climo];
   climo_wd_ptr  = new WrfData * [n_climo];
   
   
   if(!climo_prs_lvl || !climo_wd_ptr) {
      cerr << "\n\nERROR: GCPairData::set_n_climo() -> "
           << "Memory allocation error!\n\n" << flush;
      exit(1);
   }
   
   for(i=0; i<n_climo; i++) {
      climo_prs_lvl[i] = bad_data_int;
      climo_wd_ptr[i]  = (WrfData *) 0;
   }
   
   return;
}

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

void GCPairData::set_climo_prs_lvl(int i, int lvl) {

   climo_prs_lvl[i] = lvl;
   
   return;
}

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

void GCPairData::set_climo_wd_ptr(int i, WrfData *wd_ptr) {

   climo_wd_ptr[i] = wd_ptr;
   
   return;
}

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

void GCPairData::set_pd_size(int types, int masks, int interps) {
   int i, j;

   // Store the dimensions for the PairData array
   n_msg_typ = types;
   n_mask    = masks;
   n_interp  = interps;
   
   // Allocate space for the PairData array
   pd = new PairData ** [n_msg_typ];
   
   for(i=0; i<n_msg_typ; i++) {
      pd[i] = new PairData * [n_mask];
   
      for(j=0; j<n_mask; j++) {
         pd[i][j] = new PairData [n_interp];
      }
   }

   return;
}    

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

void GCPairData::set_msg_typ(int i_msg_typ, const char *name) {
   int i, j;
   
   for(i=0; i<n_mask; i++) {
      for(j=0; j<n_interp; j++) {
         pd[i_msg_typ][i][j].set_msg_typ(name);
      }
   }
   
   return;
}

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

void GCPairData::set_mask_wd(int i_mask, const char *name, WrfData *wd_ptr) {
   int i, j;
   
   for(i=0; i<n_msg_typ; i++) {
      for(j=0; j<n_interp; j++) {
         pd[i][i_mask][j].set_mask_name(name);
         pd[i][i_mask][j].set_mask_wd_ptr(wd_ptr);
      }
   }
   
   return;
}

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

void GCPairData::set_interp(int i_interp, const char *interp_mthd_str, int wdth) {
   int i, j;
   
   for(i=0; i<n_msg_typ; i++) {
      for(j=0; j<n_mask; j++) {
         pd[i][j][i_interp].set_interp_mthd(interp_mthd_str);
         pd[i][j][i_interp].set_interp_wdth(wdth);
      }
   }
   
   return;
}

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

void GCPairData::set_interp(int i_interp, int mthd, int wdth) {
   int i, j;
   
   for(i=0; i<n_msg_typ; i++) {
      for(j=0; j<n_mask; j++) {
         pd[i][j][i_interp].set_interp_mthd(mthd);
         pd[i][j][i_interp].set_interp_wdth(wdth);
      }
   }
   
   return;
}

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

void GCPairData::add_obs(float *hdr_arr, char *hdr_typ_str, 
                         char *hdr_sid_str, float *obs_arr, Grid &gr) {
   int i, j, k, x, y;
   double obs_x, obs_y, obs_prs, fcst_v, climo_v, obs_v;
   int fcst_lvl_below, fcst_lvl_above; 
   int climo_lvl_below, climo_lvl_above;

   // Check whether the grib code for the observation matches 
   // the grib code for the pair data
   if(gc.code != nint(obs_arr[3])) return;
   
   obs_prs = obs_arr[2];
   obs_v   = obs_arr[4];
      
   // Check whether the observation value contains valid data
   if(obs_v == bad_data_double) return;

   // Convert the lat/lon value to x/y
   gr.latlon_to_xy(hdr_arr[1], -1.0*hdr_arr[0], obs_x, obs_y);
   x = nint(obs_x);
   y = nint(obs_y);
   
   // Check if the observation's lat/lon is on the grid
   if(x < 0 || x >= gr.nx() || 
      y < 0 || y >= gr.ny() ) return;

   // If looking for observations on pressure levels, check whether 
   // this observation falls within the verification pressure range
   if(gc.lvl_type == PresLevel &&
      (obs_prs < gc.lvl_1 || 
       obs_prs > gc.lvl_2)) return;
   
   // If looking for observations at a vertical level, check whether 
   // the observation message type is APDSFC or SFCSHP
   if(gc.lvl_type == VertLevel &&
      strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL) return;

   // If looking for observations with an accumulation interval,
   // check whether this accumulation interval matches
   if(gc.lvl_type == AccumLevel && obs_prs != gc.lvl_1) return;

   // Find the forecast and climatology pressure levels above and below 
   // the observation point
   if(gc.lvl_type == PresLevel) {
      find_prs_lvl(1, obs_prs, fcst_lvl_below, fcst_lvl_above);
      find_prs_lvl(0, obs_prs, climo_lvl_below, climo_lvl_above);
   }
   // Set the forecast and climatology levels above and below 
   // the observation point to zero since we're looking at a single
   // level
   else if(gc.lvl_type == VertLevel || gc.lvl_type == AccumLevel) {
      fcst_lvl_below  = fcst_lvl_above  = 0;
      climo_lvl_below = climo_lvl_above = 0;
   }

   // Look through all of the PairData objects to see if the 
   // observation should be added
   
   // Check the message types
   for(i=0; i<n_msg_typ; i++) {
   
      //
      // Check for a matching PrepBufr message type
      //
      
      // Handle ANYAIR
      if(strcmp(anyair_str, pd[i][0][0].msg_typ) == 0) {
         if(strstr(anyair_msg_typ_str, hdr_typ_str) == NULL ) continue;
      }
      
      // Handle ANYSFC
      else if(strcmp(anysfc_str, pd[i][0][0].msg_typ) == 0) {
         if(strstr(anysfc_msg_typ_str, hdr_typ_str) == NULL) continue;
      }
      
      // Handle ONLYSF
      else if(strcmp(onlysf_str, pd[i][0][0].msg_typ) == 0) {
         if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL) continue;
      }
      
      // Handle all other message types
      else {
         if(strcmp(hdr_typ_str, pd[i][0][0].msg_typ) != 0) continue;
      }
   
      // Check the masking areas and points
      for(j=0; j<n_mask; j++) {
      
         // Check for the obs falling within the masking region
         if(pd[i][j][0].mask_wd_ptr != (WrfData *) 0) {
            if(!pd[i][j][0].mask_wd_ptr->s_is_on(x, y)) continue;
         }
         // Otherwise, check for the obs Station ID matching the masking SID
         else {
            if(strcmp(hdr_sid_str, pd[i][j][0].mask_name) != 0) continue;
         }
      
         // Compute the interpolated values 
         for(k=0; k<n_interp; k++) {

            // Compute the interpolated forecast value
            fcst_v = compute_interp(1, obs_x, obs_y, k,
                                    obs_prs, fcst_lvl_below, fcst_lvl_above);

            if(fcst_v == bad_data_double) continue;

            // Compute the interpolated climotological value
            climo_v = compute_interp(0, obs_x, obs_y, k, 
                                     obs_prs, climo_lvl_below, climo_lvl_above);

            // Add the forecast, climatological, and observation data
            pd[i][j][k].add_pair(fcst_v, climo_v, obs_v);

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

   return;
}

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

void GCPairData::find_prs_lvl(int fcst_flag, double prs, 
                              int &i_below, int &i_above) {
   int i, n;
   int *prs_lvl;
   double dist, dist_below, dist_above;
   
   // Check for the forecast or climo fields
   if(fcst_flag) { n = n_fcst;  prs_lvl = fcst_prs_lvl;  }
   else          { n = n_climo; prs_lvl = climo_prs_lvl; }
   
   if(n==0) {
      i_below = i_above = bad_data_int;
      return;
   }
   
   // Find the closest pressure levels above and below the observation
   dist_below = dist_above = 1.0e30;   
   for(i=0; i<n; i++) {
   
      dist = prs - prs_lvl[i];
      
      // Check for the closest level below.
      // Levels below contain higher values of pressure.
      if(dist <= 0 && abs(dist) < dist_below) { 
         dist_below = abs(dist); 
         i_below = i;
      }
      
      // Check for the closest level above.
      // Levels above contain lower values of pressure.
      if(dist >= 0 && abs(dist) < dist_above) {
         dist_above = abs(dist); 
         i_above = i;
      }
   }
   
   if(dist_below == 1.0e30 || dist_above == 1.0e30 ) {
      cerr << "\n\nERROR: GCPairData::find_prs_lvl() -> "
           << "could not find a level above and/or below the observation "
           << "pressure level of " << prs << " hp\n\n" << flush;
      exit(1);
   }

   return;
}

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

double GCPairData::compute_interp(int fcst_flag, 
                                  double obs_x, double obs_y, int i_interp,
                                  double prs, int i_below, int i_above) {
   int n;
   int *prs_lvl;
   WrfData **wd_ptr;
   double v, v_below, v_above, t;

   // Check for the forecast or climo fields
   if(fcst_flag) {
      n       = n_fcst;
      prs_lvl = fcst_prs_lvl;
      wd_ptr  = fcst_wd_ptr;
   }
   else {
      n       = n_climo;
      prs_lvl = climo_prs_lvl;
      wd_ptr  = climo_wd_ptr;
   }
   
   if(n==0) return(bad_data_double);
   
   v_below = compute_horz_interp(wd_ptr[i_below], obs_x, obs_y,
                                 pd[0][0][i_interp].interp_mthd, 
                                 pd[0][0][i_interp].interp_wdth);
   
   if(i_below == i_above) {
      v = v_below;
   }
   else {
      v_above = compute_horz_interp(wd_ptr[i_above], obs_x, obs_y,
                                    pd[0][0][i_interp].interp_mthd, 
                                    pd[0][0][i_interp].interp_wdth);
   
      // If verifying specific humidity, do vertical interpolation in 
      // the natural log of q
      if(gc.code == spfh_grib_code) {
         t = compute_vert_interp(log(v_below), prs_lvl[i_below], 
                                 log(v_above), prs_lvl[i_above], prs);
         v = exp(t);
      }
      // Vertically interpolate to the observation pressure level
      else {
         v = compute_vert_interp(v_below, prs_lvl[i_below], 
                                 v_above, prs_lvl[i_above], prs);
      }
   }
   
   return(v);
}

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

double GCPairData::compute_horz_interp(WrfData *wd_ptr, 
                                       double obs_x, double obs_y,
                                       int mthd, int wdth) {
   double v;
   int x_ll, y_ll;

   // The neighborhood width is odd, find the lower-left corner of 
   // the neighborhood
   if(wdth%2 == 1) {
      x_ll = nint(obs_x) - (wdth - 1)/2;
      y_ll = nint(obs_y) - (wdth - 1)/2;
   }
   // The neighborhood width is even, find the lower-left corner of 
   // the neighborhood
   else {
      x_ll = nint(floor(obs_x) - (wdth/2 - 1));
      y_ll = nint(floor(obs_y) - (wdth/2 - 1));
   }
   
   // Compute the interpolated value for the fields above and below
   switch(mthd) {

      case i_min:     // Minimum
         v = interp_min(wd_ptr, x_ll, y_ll, wdth, interp_threshold);
         break;
      
      case i_max:     // Maximum
         v = interp_max(wd_ptr, x_ll, y_ll, wdth, interp_threshold);
         break;
      
      case i_median:  // Median
         v = interp_median(wd_ptr, x_ll, y_ll, wdth, interp_threshold);
         break;
      
      case i_uw_mean: // Unweighted Mean
         v = interp_uw_mean(wd_ptr, x_ll, y_ll, wdth, interp_threshold);
         break;
      
      case i_dw_mean: // Distance-Weighted Mean
         v = interp_dw_mean(wd_ptr, x_ll, y_ll, wdth, obs_x, obs_y, 
                            dw_mean_pow, interp_threshold);
         break;

      case i_ls_fit:  // Least-squares fit
         v = interp_ls_fit(wd_ptr, x_ll, y_ll, wdth, obs_x, obs_y, 
                           interp_threshold);
         break;
 
      default:
         cerr << "\n\nERROR: GCPairData::compute_horz_interp() -> "
              << "unexpected interpolation method encountered: "
              << mthd << "\n\n" << flush;
         exit(1);
         break;
   }
   
   return(v);
}

////////////////////////////////////////////////////////////////////////
//
// Vertically interpolate between values "v1" and "v2" at pressure levels
// "prs1" and "prs2" to pressure level "to_prs".
//
////////////////////////////////////////////////////////////////////////

double GCPairData::compute_vert_interp(double v1, double prs1, 
                                       double v2, double prs2,
                                       double to_prs) {
   double v_interp;

   if(prs1 <= 0.0 || prs2 <= 0.0 || to_prs <= 0.0) {
      cerr << "\n\nERROR: GCPairData::compute_vert_interp() -> "
           << "pressure shouldn't be <= zero!\n\n" << flush;
      exit(1);
   }
   
   v_interp = v1 + ((v2-v1)*log(prs1/to_prs)/log(prs1/prs2));
   
   return(v_interp);
}

////////////////////////////////////////////////////////////////////////
//
// Code for class CIInfo
//
////////////////////////////////////////////////////////////////////////

CIInfo::CIInfo() {
   init_from_scratch();
}

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

CIInfo::~CIInfo() {
   clear();
}

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

CIInfo::CIInfo(const CIInfo &c) {
   
   init_from_scratch();

   assign(c);
}

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

CIInfo & CIInfo::operator=(const CIInfo &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

void CIInfo::init_from_scratch() {

   v_cl = (double *) 0;
   v_cu = (double *) 0;

   clear();

   return;
}

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

void CIInfo::clear() {
   
   n = 0;
   v = 0.0;
   if(v_cl) { delete [] v_cl; v_cl = (double *) 0; }
   if(v_cu) { delete [] v_cu; v_cu = (double *) 0; }

   return;
}

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

void CIInfo::assign(const CIInfo &c) {
   int i;

   clear();
   
   allocate_ci(c.n);
 
   v = c.v;
   
   for(i=0; i<c.n; i++) {
      v_cl[i] = c.v_cl[i];
      v_cu[i] = c.v_cu[i];
   }

   return;
}

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

void CIInfo::allocate_ci(int i) {
   
   n = i;
   
   if(n > 0) {
      v_cl = new double [n];
      v_cu = new double [n];

      if(!v_cl || !v_cu) {
         cerr << "\n\nERROR: CIInfo::allocate_ci() -> "
              << "Memory allocation error!\n\n" << flush;
         exit(1);
      }
   }

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class CTSInfo
//
////////////////////////////////////////////////////////////////////////

CTSInfo::CTSInfo() {
   init_from_scratch();
}

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

CTSInfo::~CTSInfo() {
   clear();
}

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

CTSInfo::CTSInfo(const CTSInfo &c) {
   
   init_from_scratch();

   assign(c);
}

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

CTSInfo & CTSInfo::operator=(const CTSInfo &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

void CTSInfo::init_from_scratch() {

   name  = (char *) 0;
   alpha = (double *) 0;

   clear();

   return;
}

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

void CTSInfo::clear() {
   
   cts.zero_out();
   cts_thresh = 0.0;
   cts_thresh_ind = 0;
   
   n_alpha = 0;
   if(alpha) { delete [] alpha; alpha = (double *) 0; }
   if(name)  { delete [] name;  name = (char *) 0; }
   
   baser.clear();
   fmean.clear();
   acc.clear();
   fbias.clear();
   pody.clear();
   podn.clear();
   pofd.clear();
   far.clear();
   csi.clear();
   gss.clear();
   hk.clear();
   hss.clear();
   odds.clear();

   return;
}

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

void CTSInfo::assign(const CTSInfo &c) {
   int i;
   
   clear();
   
   cts = c.cts;
   cts_thresh = c.cts_thresh;
   cts_thresh_ind = c.cts_thresh_ind;
   
   allocate_n_alpha(c.n_alpha);
   for(i=0; i<c.n_alpha; i++) { alpha[i] = c.alpha[i]; }
   
   set_name(c.name);
   
   baser = c.baser;
   fmean = c.fmean;
   acc = c.acc;
   fbias = c.fbias;
   pody = c.pody;
   podn = c.podn;
   pofd = c.pofd;
   far = c.far;
   csi = c.csi;
   gss = c.gss;
   hk = c.hk;
   hss = c.hss;
   odds = c.odds;

   return;
}

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

void CTSInfo::allocate_n_alpha(int i) {

   n_alpha = i;

   if(n_alpha > 0) {

      alpha = new double [n_alpha];

      if(!alpha) {
         cerr << "\n\nERROR: CTSInfo::allocate_n() -> "
              << "Memory allocation error!\n\n" << flush;
         exit(1);
      }
   
      baser.allocate_ci(n_alpha);
      fmean.allocate_ci(n_alpha);
      acc.allocate_ci(n_alpha);
      fbias.allocate_ci(n_alpha);
      pody.allocate_ci(n_alpha);
      podn.allocate_ci(n_alpha);
      pofd.allocate_ci(n_alpha);
      far.allocate_ci(n_alpha);
      csi.allocate_ci(n_alpha);
      gss.allocate_ci(n_alpha);
      hk.allocate_ci(n_alpha);
      hss.allocate_ci(n_alpha);
      odds.allocate_ci(n_alpha);
   }

   return;
}

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

void CTSInfo::set_name(const char *c) {
   
   if(c) {
      name = new char [strlen(c)+1];
      strcpy(name, c);
   }
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class CNTInfo
//
////////////////////////////////////////////////////////////////////////

CNTInfo::CNTInfo() {
   init_from_scratch();
}

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

CNTInfo::~CNTInfo() {
   clear();
}

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

CNTInfo::CNTInfo(const CNTInfo &c) {
   
   init_from_scratch();

   assign(c);
}

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

CNTInfo & CNTInfo::operator=(const CNTInfo &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

void CNTInfo::init_from_scratch() {

   name  = (char *) 0;
   alpha = (double *) 0;

   clear();

   return;
}

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

void CNTInfo::clear() {

   n = 0;
   n_alpha = 0;   

   if(name)  { delete [] name;   name = (char *) 0; }   
   if(alpha) { delete [] alpha; alpha = (double *) 0; }
   
   fbar.clear();
   fstdev.clear();
   obar.clear();
   ostdev.clear();
   pr_corr.clear();
   me.clear();
   estdev.clear();

   sp_corr = kt_corr = 0.0;

   n_ranks = frank_ties = orank_ties = 0;
   bias = mbias = mae = mse = bcmse = rmse = 0.0;
   e10 = e25 = e50 = e75 = e90 = 0.0;

   return;
}

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

void CNTInfo::assign(const CNTInfo &c) {
   int i;

   clear();
   
   n = c.n;
   allocate_n_alpha(c.n_alpha);
   for(i=0; i<c.n; i++) { alpha[i] = c.alpha[i]; }
      
   set_name(c.name);

   fbar        = c.fbar;
   fstdev      = c.fstdev;
   obar        = c.obar;
   ostdev      = c.ostdev;
   pr_corr     = c.pr_corr;
   sp_corr     = c.sp_corr;
   kt_corr     = c.kt_corr;
   n_ranks     = c.n_ranks;
   frank_ties = c.frank_ties;
   orank_ties  = c.orank_ties;
   me          = c.me;
   estdev      = c.estdev;
   
   bias        = c.bias;
   mbias       = c.mbias;
   mae         = c.mae;
   mse         = c.mse;
   bcmse       = c.bcmse;
   rmse        = c.rmse;
   e10         = c.e10;
   e25         = c.e25;
   e50         = c.e50;
   e75         = c.e75;
   e90         = c.e90;

   return;
}

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

void CNTInfo::allocate_n_alpha(int i) {

   n_alpha = i;

   if(n_alpha > 0) {

      alpha = new double [n_alpha];

      if(!alpha) {
         cerr << "\n\nERROR: CNTInfo::allocate_n_alpha() -> "
              << "Memory allocation error!\n\n" << flush;
         exit(1);
      }
   
      fbar.allocate_ci(n_alpha);
      fstdev.allocate_ci(n_alpha);
      obar.allocate_ci(n_alpha);
      ostdev.allocate_ci(n_alpha);
      pr_corr.allocate_ci(n_alpha);
      me.allocate_ci(n_alpha);
      estdev.allocate_ci(n_alpha);
   }

   return;
}

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

void CNTInfo::set_name(const char *c) {
   
   if(c) {
      name = new char [strlen(c)+1];
      strcpy(name, c);
   }
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class SL1L2Info
//
////////////////////////////////////////////////////////////////////////

SL1L2Info::SL1L2Info() {
   init_from_scratch();
}

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

SL1L2Info::~SL1L2Info() {
   clear();
}

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

SL1L2Info::SL1L2Info(const SL1L2Info &c) {
   
   init_from_scratch();

   assign(c);
}

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

SL1L2Info & SL1L2Info::operator=(const SL1L2Info &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

SL1L2Info & SL1L2Info::operator+=(const SL1L2Info &c) {
   SL1L2Info s_info;
   
   s_info.scount  = scount  + c.scount;
   
   if(s_info.scount > 0) {
      s_info.fbar  = (fbar*scount  + c.fbar*c.scount) /s_info.scount;
      s_info.obar  = (obar*scount  + c.obar*c.scount) /s_info.scount;
      s_info.fobar = (fobar*scount + c.fobar*c.scount)/s_info.scount;
      s_info.ffbar = (ffbar*scount + c.ffbar*c.scount)/s_info.scount;
      s_info.oobar = (oobar*scount + c.oobar*c.scount)/s_info.scount;
   }
   
   s_info.sacount  = sacount  + c.sacount;
   
   if(s_info.sacount > 0) {
      s_info.fabar  = (fabar*sacount  + c.fabar*c.sacount) /s_info.sacount;
      s_info.oabar  = (oabar*sacount  + c.oabar*c.sacount) /s_info.sacount;
      s_info.foabar = (foabar*sacount + c.foabar*c.sacount)/s_info.sacount;
      s_info.ffabar = (ffabar*sacount + c.ffabar*c.sacount)/s_info.sacount;
      s_info.ooabar = (ooabar*sacount + c.ooabar*c.sacount)/s_info.sacount;
   }
   
   assign(s_info);

   return(*this);
}

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

void SL1L2Info::init_from_scratch() {

   clear();

   return;
}

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

void SL1L2Info::clear() {

   // SL1L2 Quantities
   fbar   = obar   = 0.0;
   fobar  = ffbar  = oobar  = 0.0;
   scount = 0;
   
   // SAL1L2 Quantities
   fabar  = oabar  = 0.0;
   foabar = ffabar = ooabar = 0.0;
   sacount = 0;

   return;
}

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

void SL1L2Info::assign(const SL1L2Info &c) {

   clear();
   
   // SL1L2 Quantities
   fbar    = c.fbar;
   obar    = c.obar;
   fobar   = c.fobar;
   ffbar   = c.ffbar;
   oobar   = c.oobar;
   scount  = c.scount;
   
   // SAL1L2 Quantities
   fabar   = c.fabar;
   oabar   = c.oabar;
   foabar  = c.foabar;
   ffabar  = c.ffabar;
   ooabar  = c.ooabar;
   sacount = c.sacount;
   
   return;
}

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

void compute_cnt(SL1L2Info &s, int aflag, double alpha, CNTInfo &c) {
   double den, f, o;
   double cv_normal_l, cv_normal_u, cv_chi2_l, cv_chi2_u;
   double v, cl, cu;

   // Handle the count
   if(!aflag) c.n = s.scount;
   else       c.n = s.sacount;

   // Range check the alpha value supplied
   if(alpha <= 0 || alpha >= 1.0) {
      cerr << "\n\nERROR: compute_cnt() -> " 
           << "Confidence interval alpha values (" 
           <<  alpha << ") must be greater than 0 "
           << "and less than 1.\n\n" << flush;
      exit(1);
   }
   
   // Allocate enough space and store all of the alpha values to be used
   c.n_alpha = 1;
   c.allocate_n_alpha(c.n_alpha);
   c.alpha[0] = alpha;
   
   // Set the quantities that can't be derived from SL1L2Info to bad data
   c.sp_corr    = bad_data_double;
   c.kt_corr    = bad_data_double;
   c.n_ranks    = bad_data_int;
   c.frank_ties = bad_data_int;
   c.orank_ties = bad_data_int;
   c.mae        = bad_data_double;
   c.e10        = bad_data_double;
   c.e25        = bad_data_double;
   c.e50        = bad_data_double;
   c.e75        = bad_data_double;
   c.e90        = bad_data_double;
   
   // Compute the forecast mean
   if(!aflag) c.fbar.v = s.fbar;
   else       c.fbar.v = s.fabar;

   // Compute the forecast standard deviation
   if(!aflag) c.fstdev.v = compute_stdev(s.fbar*s.scount, s.ffbar*s.scount, s.scount);
   else       c.fstdev.v = compute_stdev(s.fabar*s.sacount, s.ffabar*s.sacount, s.sacount);

   // Compute the observation mean
   if(!aflag) c.obar.v = s.obar;
   else       c.obar.v = s.oabar;

   // Compute the observation standard deviation
   if(!aflag) c.ostdev.v = compute_stdev(s.obar*s.scount, s.fobar*s.scount, s.scount);
   else       c.ostdev.v = compute_stdev(s.oabar*s.sacount, s.foabar*s.sacount, s.sacount);

   // Compute the f*o mean
   if(!aflag) c.fobar = s.fobar;
   else       c.fobar = s.foabar;

   // Compute the forecast squared mean
   if(!aflag) c.ffbar = s.ffbar;
   else       c.ffbar = s.ffabar;

   // Compute the observation squared mean
   if(!aflag) c.oobar = s.oobar;
   else       c.oobar = s.ooabar;

   // Compute the bias
   c.bias = c.fbar.v - c.obar.v;
   
   // Compute the multiplicative bias
   if(c.obar.v == 0.0) c.mbias = bad_data_double;
   else                c.mbias = c.fbar.v/c.obar.v;
   
   // Compute the correlation coefficient
   den = sqrt((c.n*c.ffbar*c.n - c.fbar.v*c.n*c.fbar.v*c.n)*
              (c.n*c.oobar*c.n - c.obar.v*c.n*c.obar.v*c.n));
   if(den == 0.0) c.pr_corr.v = bad_data_double; 
   else           c.pr_corr.v = ((c.n*c.fobar*c.n) - (c.fbar.v*c.n*c.obar.v*c.n))/den;

   // Compute mean error
   c.me.v = c.fbar.v - c.obar.v;

   // Compute the mean squared error
   c.mse = c.ffbar + c.oobar - 2.0*c.fobar;

   // Compute the standard deviation of the mean error
   c.estdev.v = compute_stdev(c.me.v*c.n, c.mse*c.n, c.n);
   
   // Compute the bias corrected mean squared error (decomposition of MSE)
   f = c.fbar.v;
   o = c.obar.v;
   c.bcmse = c.mse - (f-o)*(f-o);

   // Compute root mean squared error
   c.rmse = sqrt(c.mse);

   //
   // Compute confidence intervals for the alpha value specified
   // 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.
   //
      
   // Check for the degenerate case
   if(c.n <= 1) {
      c.fbar.v_cl[0]    = c.fbar.v_cu[0]    = bad_data_double;
      c.fstdev.v_cl[0]  = c.fstdev.v_cu[0]  = bad_data_double;
      c.obar.v_cl[0]    = c.obar.v_cu[0]    = bad_data_double;
      c.ostdev.v_cl[0]  = c.ostdev.v_cu[0]  = bad_data_double;
      c.pr_corr.v_cl[0] = c.pr_corr.v_cu[0] = bad_data_double;
      c.me.v_cl[0]      = c.me.v_cu[0]      = bad_data_double;
      c.estdev.v_cl[0]  = c.estdev.v_cu[0]  = bad_data_double;
      return;
   }
   
   // Compute the critical values for the Normal or Student's-T distribution
   // based on the sample size
   if(c.n >= large_sample_threshold) {
      cv_normal_l = normal_cdf_inv(c.alpha[0]/2.0, 0.0, 1.0);
      cv_normal_u = normal_cdf_inv(1.0 - (c.alpha[0]/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(c.alpha[0]/2.0, c.n - 1);
      cv_normal_u = students_t_cdf_inv(1.0 - (c.alpha[0]/2.0), c.n - 1);
   }

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

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

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

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

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

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

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

   // Compute confidence interval for the error standard deviation
   c.estdev.v_cl[0] = sqrt((c.n - 1)*c.estdev.v*c.estdev.v/cv_chi2_u);
   c.estdev.v_cu[0] = sqrt((c.n - 1)*c.estdev.v*c.estdev.v/cv_chi2_l);
   
   return;
}

////////////////////////////////////////////////////////////////////////
//
// Code for class VL1L2Info
//
////////////////////////////////////////////////////////////////////////

VL1L2Info::VL1L2Info() {
   init_from_scratch();
}

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

VL1L2Info::~VL1L2Info() {
   clear();
}

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

VL1L2Info::VL1L2Info(const VL1L2Info &c) {
   
   init_from_scratch();

   assign(c);
}

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

VL1L2Info & VL1L2Info::operator=(const VL1L2Info &c) {

   if(this == &c) return(*this);

   assign(c);

   return(*this);
}

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

VL1L2Info & VL1L2Info::operator+=(const VL1L2Info &c) {
   VL1L2Info v_info;

   v_info.vcount  = vcount  + c.vcount;
   
   if(v_info.vcount > 0) {
      v_info.ufbar   = (ufbar*vcount   + c.ufbar*c.vcount)  /v_info.vcount;
      v_info.vfbar   = (vfbar*vcount   + c.vfbar*c.vcount)  /v_info.vcount;
      v_info.uobar   = (uobar*vcount   + c.uobar*c.vcount)  /v_info.vcount;
      v_info.vobar   = (vobar*vcount   + c.vobar*c.vcount)  /v_info.vcount;
      v_info.uvfobar = (uvfobar*vcount + c.uvfobar*c.vcount)/v_info.vcount;
      v_info.uvffbar = (uvffbar*vcount + c.uvffbar*c.vcount)/v_info.vcount;
      v_info.uvoobar = (uvoobar*vcount + c.uvoobar*c.vcount)/v_info.vcount;
   }
   
   v_info.vacount  = vacount  + c.vacount;
   
   if(v_info.vacount > 0) {
      v_info.ufabar   = (ufabar*vacount   + c.ufabar*c.vacount)  /v_info.vacount;
      v_info.vfabar   = (vfabar*vacount   + c.vfabar*c.vacount)  /v_info.vacount;
      v_info.uoabar   = (uoabar*vacount   + c.uoabar*c.vacount)  /v_info.vacount;
      v_info.voabar   = (voabar*vacount   + c.voabar*c.vacount)  /v_info.vacount;
      v_info.uvfoabar = (uvfoabar*vacount + c.uvfoabar*c.vacount)/v_info.vacount;
      v_info.uvffabar = (uvffabar*vacount + c.uvffabar*c.vacount)/v_info.vacount;
      v_info.uvooabar = (uvooabar*vacount + c.uvooabar*c.vacount)/v_info.vacount;
   }
      
   assign(v_info);

   return(*this);
}

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

void VL1L2Info::init_from_scratch() {

   clear();

   return;
}

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

void VL1L2Info::clear() {

   // VL1L2 Quantities
   ufbar    = vfbar    = uobar    = vobar   = 0.0;
   uvfobar  = uvffbar  = uvoobar  = 0.0;
   vcount   = 0;
   
   // VAL1L2 Quantities
   ufabar   = vfabar   = uoabar   = voabar  = 0.0;
   uvfoabar = uvffabar = uvooabar = 0.0;
   vacount  = 0;

   return;
}

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

void VL1L2Info::assign(const VL1L2Info &c) {

   clear();
   
   // VL1L2 Quantities
   ufbar    = c.ufbar;
   vfbar    = c.vfbar;
   uobar    = c.uobar;
   vobar    = c.vobar;
   uvfobar  = c.uvfobar;
   uvffbar  = c.uvffbar;
   uvoobar  = c.uvoobar;
   vcount   = c.vcount;
   
   // VAL1L2 Quantities
   ufabar   = c.ufabar;
   vfabar   = c.vfabar;
   uoabar   = c.uoabar;
   voabar   = c.voabar;
   uvfoabar = c.uvfoabar;
   uvffabar = c.uvffabar;
   uvooabar = c.uvooabar;
   vacount  = c.vacount;

   return;
}

////////////////////////////////////////////////////////////////////////
//
// Begin code for misc functions
//
////////////////////////////////////////////////////////////////////////

void parse_gcinfo(const char *gc_str, GCInfo &gc) {
   char tmp_str[max_str_len], tmp2_str[max_str_len];
   char *ptr, *ptr2;
   int j;
      
   // Initialize the temp string
   strcpy(tmp_str, gc_str);
   
   ptr = strchr(tmp_str, '/');
   if(ptr == NULL) {
      cerr << "\n\nERROR: parse_gcinfo() -> "
           << "each grib code specified must be followed by an "
           << "accumulation, level, or presssure level indicator \"" 
           << gc_str << "\".\n\n" << flush;
      exit(1);
   }
   
   // Copy the grib code into another string and null terminate it
   for(j=0; tmp_str[j] != '/'; j++) tmp2_str[j] = tmp_str[j];
   tmp2_str[j] = '\0';
   gc.code = str_to_grib_code(tmp2_str);
   
   // Advance the pointer past the '/'
   ptr++;
   if(*ptr != 'A' && *ptr != 'Z' && *ptr != 'P') {
      cerr << "\n\nERROR: parse_gcinfo() -> "
           << "each grib code specified must be followed by an "
           << "accumulation, level, or presssure level indicator \"" 
           << gc_str << "\".\n\n" << flush;
      exit(1);   
   }
   
   // Set the level type
   if(      *ptr == 'A') gc.lvl_type = AccumLevel;
   else if (*ptr == 'Z') gc.lvl_type = VertLevel;
   else if (*ptr == 'P') gc.lvl_type = PresLevel;
   else                  gc.lvl_type = NoLevel;
   
   // Store the level string
   gc.set_lvl_str(ptr);
      
   // Advance the pointer past the 'A', 'Z' or 'P'
   ptr++;
   gc.lvl_1 = atoi(ptr);
   
   // Look for a '-' and a second level indicator
   ptr2 = strchr(ptr, '-');
   if(ptr2 != NULL) {
      gc.lvl_2 = atoi(++ptr2);
   }
   else {
      gc.lvl_2 = gc.lvl_1;
   }
   
   if(gc.lvl_type != PresLevel &&
      gc.lvl_1    != gc.lvl_2) {
      cerr << "\n\nERROR: parse_gcinfo() -> "
           << "ranges of values are only supported for pressure levels \""
           << gc_str << "\".\n\n" << flush;
      exit(1);
   }
   
   if(gc.lvl_1 > gc.lvl_2) {
      cerr << "\n\nERROR: parse_gcinfo() -> "
           << "the first level must be less than the second level \""
           << gc_str << "\".\n\n" << flush;
      exit(1);
   }   
   
   return;
}

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

int parse_thresh(const char *thresh_str, int *&n_thresh_arr, 
                 double **&thresh_arr, int **&thresh_ind_arr) {
   char tmp_str[max_str_len], **temp_sub_str;
   char *c = (char *) 0;
   int i, n;
   
   // Compute the number of tokens in the string based on "|"
   n = num_tokens(thresh_str, "|");
   
   // Check for no tokens in string
   if(n == 0) return(0);
   
   // Allocate space
   temp_sub_str = new char * [n];
   n_thresh_arr = new int [n];
   thresh_arr = new double * [n];
   thresh_ind_arr = new int * [n];
   for(i=0; i<n; i++) thresh_arr[i] = (double *) 0;
   for(i=0; i<n; i++) temp_sub_str[i] = new char [max_str_len];
      
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, thresh_str);

   // Tokenize the string
   c = strtok(tmp_str, "|");
   strcpy(temp_sub_str[0], c);
   
   // Parse remaining tokens
   for(i=1; i<n; i++) strcpy(temp_sub_str[i], strtok(0, "|"));
      
   // Parse the sub strings
   for(i=0; i<n; i++) {
      n_thresh_arr[i] = parse_thresh_entry(temp_sub_str[i], 
                                           thresh_arr[i],
                                           thresh_ind_arr[i]);
      if(n_thresh_arr[i] <= 0) {
         cerr << "\n\nERROR: parse_thresh() -> "
              << "at least one threshold value must be provided for each grib code "
              << "listed in vx_grib_code. \"" << thresh_str << "\".\n\n" << flush;
         exit(1);      
      }
   }
 
   for(i=0; i<n; i++) { 
      if(temp_sub_str[i]) { 
         delete [] (temp_sub_str[i]);
         temp_sub_str[i] = (char *) 0;
      }
   }
   if(temp_sub_str) {
      delete [] (temp_sub_str);
      temp_sub_str = (char **) 0;
   }
   
   return(n);
}

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

int parse_thresh_entry(const char *dbl_str, double *&dbl_arr, int *&dbl_ind_arr) {
   char tmp_str[max_str_len], a, b;
   char *c = (char *) 0;
   int n, i;

   // Compute the number of tokens in the string based on " "
   n = num_tokens(dbl_str, " ");
   
   // Check for no tokens in string
   if(n == 0) return(0);
   
   // Allocate space for the list of tokens
   dbl_arr = new double [n];
   dbl_ind_arr = new int [n];
   
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, dbl_str);

   // Tokenize the string and store the threshold values
   c = strtok(tmp_str, " ");
   
   // Check that the threshold string is more than 2 characters long
   if(strlen(c) < 3) {
      cerr << "\n\nERROR: parse_thresh_entry() -> "
           << "each threshold value must be preceeded by one of "
           << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
           << "thresholding to be performed \"" << dbl_str
           << "\".\n\n" << flush;
      exit(1);
   }
   
   a = tolower(c[0]);
   b = tolower(c[1]);
   
   if     (a == 'l' && b == 't') dbl_ind_arr[0] = thresh_lt; // less than
   else if(a == 'l' && b == 'e') dbl_ind_arr[0] = thresh_le; // less than or equal to
   else if(a == 'e' && b == 'q') dbl_ind_arr[0] = thresh_eq; // equal to
   else if(a == 'n' && b == 'e') dbl_ind_arr[0] = thresh_ne; // not equal to
   else if(a == 'g' && b == 't') dbl_ind_arr[0] = thresh_gt; // greater than
   else if(a == 'g' && b == 'e') dbl_ind_arr[0] = thresh_ge; // greater than or equal to
   else {
      cerr << "\n\nERROR: parse_thresh_entry() -> "
           << "each threshold value must be preceeded by one of "
           << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
           << "thresholding to be performed \"" << dbl_str
           << "\".\n\n" << flush;
      exit(1);   
   }
   c += 2;
   dbl_arr[0] = atof(c);
   
   // Parse remaining tokens
   for(i=1; i<n; i++) {
      
      c = strtok(0, " ");
      
      // Check that the threshold string is more than 2 characters long
      if(strlen(c) < 3) {
         cerr << "\n\nERROR: parse_thresh_entry() -> "
              << "each threshold value must be preceeded by one of "
              << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
              << "thresholding to be performed \"" << dbl_str
              << "\".\n\n" << flush;
         exit(1);
      }
      
      a = tolower(c[0]);
      b = tolower(c[1]);
   
      if     (a == 'l' && b == 't') dbl_ind_arr[i] = thresh_lt; // less than
      else if(a == 'l' && b == 'e') dbl_ind_arr[i] = thresh_le; // less than or equal to
      else if(a == 'e' && b == 'q') dbl_ind_arr[i] = thresh_eq; // equal to
      else if(a == 'n' && b == 'e') dbl_ind_arr[i] = thresh_ne; // not equal to
      else if(a == 'g' && b == 't') dbl_ind_arr[i] = thresh_gt; // greater than
      else if(a == 'g' && b == 'e') dbl_ind_arr[i] = thresh_ge; // greater than or equal to
      else {
         cerr << "\n\nERROR: parse_thresh_entry() -> "
              << "each threshold value must be preceeded by one of "
              << "\"lt, le, eq, ne, gt, ge\" to indicate the type of "
              << "thresholding to be performed \"" << dbl_str
              << "\".\n\n" << flush;
         exit(1);   
      }
      c += 2;
      dbl_arr[i] = atof(c);
   }
   
   return(n);  
}

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

int parse_message_type(const char *msg_typ_str, char **&msg_typ_arr) {
   char tmp_str[max_str_len];
   char *c = (char *) 0;
   int n, i;

   // Compute the number of tokens in the string based on " "
   n = num_tokens(msg_typ_str, " ");
   
   // Check for no tokens in string
   if(n == 0) return(0);
   
   // Allocate space for the list of tokens
   msg_typ_arr = new char * [n];
   
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, msg_typ_str);

   // Tokenize the string and store the double values
   c = strtok(tmp_str, " ");
   msg_typ_arr[0] = new char [strlen(c)+1];
   strcpy(msg_typ_arr[0], c);
   
   // Parse remaining tokens
   for(i=1; i<n; i++) {
      c = strtok(0, " ");
      msg_typ_arr[i] = new char [strlen(c)+1];
      strcpy(msg_typ_arr[i], c);
   }
   
   return(n);
}

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

int parse_dbl_list(const char *dbl_str, double *&dbl_arr) {
   char tmp_str[max_str_len];
   char *c = (char *) 0;
   int n, i;

   // Compute the number of tokens in the string based on " "
   n = num_tokens(dbl_str, " ");
   
   // Check for no tokens in string
   if(n == 0) return(0);
   
   // Allocate space for the list of tokens
   dbl_arr = new double [n];   
   
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, dbl_str);

   // Tokenize the string and store the double values
   c = strtok(tmp_str, " ");
   dbl_arr[0] = atof(c);
   
   // Parse remaining tokens
   for(i=1; i<n; i++) dbl_arr[i] = atof(strtok(0, " "));
   
   return(n);  
}

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

int parse_int_list(const char *int_str, int *&int_arr) {
   char tmp_str[max_str_len];
   char *c = (char *) 0;
   int n, i;

   // Compute the number of tokens in the string based on " "
   n = num_tokens(int_str, " ");
   
   // Check for no tokens in string
   if(n == 0) return(0);
   
   // Allocate space for the list of tokens
   int_arr = new int [n];   
   
   // Initialize the temp string for use in tokenizing
   strcpy(tmp_str, int_str);

   // Tokenize the string and store the integer values
   c = strtok(tmp_str, " ");
   int_arr[0] = nint(atof(c));
   
   // Parse remaining tokens
   for(i=1; i<n; i++) int_arr[i] = nint(atof(strtok(0, " ")));
   
   return(n);  
}

////////////////////////////////////////////////////////////////////////
            
int max_int(const int *v_int, int n) {
   int i, v_max;

   if(n <= 0) return(0);
   
   v_max = v_int[0];
   for(i=1; i<n; i++) if(v_int[i] > v_max) v_max = v_int[i];

   return(v_max);
}

////////////////////////////////////////////////////////////////////////
            
int min_int(const int *v_int, int n) {
   int i, v_min;

   if(n <= 0) return(0);
   
   v_min = v_int[0];
   for(i=1; i<n; i++) if(v_int[i] < v_min) v_min = v_int[i];

   return(v_min);
}

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

void dbl_to_str(double v, char *v_str) {

   dbl_to_str(v, v_str, default_precision);

   return; 
}

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

void dbl_to_str(double v, char *v_str, int precision) {
   char fmt_str[32];
   
   sprintf(fmt_str, "%s%i%s", "%.", precision, "f");

   if(v == bad_data_double) sprintf(v_str, "%i", bad_data_int);
   else                     sprintf(v_str, fmt_str, v);

   return; 
}

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

double compute_stdev(double sum, double sum_sq, int n) {
   double sigma;
   
   if(n <= 1) {
      sigma = bad_data_double;
   }
   else {
      sigma = sqrt((sum_sq - sum*sum/(double) n)/((double) (n - 1)));
   }
   
   return(sigma);
}

////////////////////////////////////////////////////////////////////////
//
// For the following format and header string commands, the character
// pointer (ptr) must contain enough allocated memory to store the 
// resulting string.
//
////////////////////////////////////////////////////////////////////////

void fho_fmt_str(char *fmt_str) {

   // Setup the FHO format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void fho_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the FHO format string
   fho_fmt_str(fmt_str);    

   // Setup the FHO header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "FHO_T", "VAR", "LEVEL", 
           "TOTAL", "F_RATE", "H_RATE", "O_RATE", 
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void ctc_fmt_str(char *fmt_str) {

   // Setup the CTC format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void ctc_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the CTC format string
   ctc_fmt_str(fmt_str);    

   // Setup the CTC header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "CTC_T", "VAR", "LEVEL", 
           "TOTAL", "FY_OY", "FY_ON", "FN_OY", "FN_ON", 
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void ctp_fmt_str(char *fmt_str) {

   // Setup the CTP format string
   sprintf(fmt_str, "%s %s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s");

   return;
}

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

void ctp_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the CTP format string
   ctp_fmt_str(fmt_str);    

   // Setup the CTP header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "CTP_T", "VAR", "LEVEL", 
           "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");

   return;
}

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

void cfp_fmt_str(char *fmt_str) {

   // Setup the CFP format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void cfp_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the CFP format string
   cfp_fmt_str(fmt_str);    

   // Setup the CFP header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "CFP_T", "VAR", "LEVEL", 
           "TOTAL", "FY_OY_FP", "FY_ON_FP", "FN_OY_FP", "FN_ON_FP", 
           "FY", "FN",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void cop_fmt_str(char *fmt_str) {

   // Setup the COP format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void cop_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the COP format string
   cop_fmt_str(fmt_str);    

   // Setup the COP header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "COP_T", "VAR", "LEVEL", 
           "TOTAL", "FY_OY_OP", "FY_ON_OP", "FN_OY_OP", "FN_ON_OP", 
           "OY", "ON",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void cts_fmt_str(char *fmt_str) {

   // Setup the CTS format string
   sprintf(fmt_str, "%s %s %s %s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void cts_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the CTS format string
   cts_fmt_str(fmt_str);    

   // Setup the CTS header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "CTS_T/ALPHA", "VAR", "LEVEL", 
           "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");

   return;
}

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

void cnt_fmt_str(char *fmt_str) {

   // Setup the CNT format string
   sprintf(fmt_str, "%s %s %s %s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s",
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void cnt_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the CNT format string
   cnt_fmt_str(fmt_str);    

   // Setup the CNT header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "CNT/ALPHA", "VAR", "LEVEL", 
           "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");

   return;
}

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

void sl1l2_fmt_str(char *fmt_str) {

   // Setup the SL1L2 format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void sl1l2_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the SL1L2 format string
   sl1l2_fmt_str(fmt_str);    

   // Setup the SL1L2 header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "SL1L2", "VAR", "LEVEL", 
           "TOTAL",
           "FBAR", "OBAR", "FOBAR", "FFBAR", "OOBAR",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void sal1l2_fmt_str(char *fmt_str) {

   // Setup the SAL1L2 format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void sal1l2_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the SAL1L2 format string
   sal1l2_fmt_str(fmt_str);    

   // Setup the SAL1L2 header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "SAL1L2", "VAR", "LEVEL", 
           "TOTAL",
           "FABAR", "OABAR", "FOABAR", "FFABAR", "OOABAR",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void vl1l2_fmt_str(char *fmt_str) {

   // Setup the VL1L2 format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void vl1l2_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the VL1L2 format string
   vl1l2_fmt_str(fmt_str);    

   // Setup the VL1L2 header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "VL1L2", "VAR", "LEVEL", 
           "TOTAL",
           "UFBAR", "VFBAR", "UOBAR", "VOBAR", 
           "UVFOBAR", "UVFFBAR", "UVOOBAR",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

void val1l2_fmt_str(char *fmt_str) {

   // Setup the VAL1L2 format string
   sprintf(fmt_str, "%s %s", 
           hdr_fmt_str,
           "%-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s %-15s");

   return;
}

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

void val1l2_hdr_str(char *hdr_str) {
   char fmt_str[max_str_len];
   
   // Retrieve the VAL1L2 format string
   val1l2_fmt_str(fmt_str);    

   // Setup the VAL1L2 header line  
   sprintf(hdr_str, fmt_str, 
           "VRS", "MODEL", "FCST_LEAD", "FCST_VALID", "OBTYPE", 
           "VX_MASK", "VAL1L2", "VAR", "LEVEL", 
           "TOTAL",
           "UFABAR", "VFABAR", "UOABAR", "VOABAR", 
           "UVFOABAR", "UVFFABAR", "UVOOABAR",
           "INTERP_MTHD", "INTERP_PNTS");

   return;
}

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

double interp_min(WrfData *wd, int x_ll, int y_ll, int wdth, double t) {
   int x, y, count;
   double v, min_v;
   
   // Search the neighborhood for the minimum value
   count = 0;
   min_v = 1.0e30;
   for(x=x_ll; x<x_ll+wdth; x++) {
      if(x < 0 || x >= wd->get_nx()) continue;
   
      for(y=y_ll; y<y_ll+wdth; y++) {
         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         v = wd->get_xy_double(x, y);
         if(v < min_v) min_v = v;
         count++;
      } // end for y
   } // end for x
   
   // Check whether enough valid grid points were found to trust
   // the minimum value computed
   if( (double) count/(wdth*wdth) < t
    || count == 0) {
      min_v = bad_data_double;
   }
   
   return(min_v);
}

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

double interp_max(WrfData *wd, int x_ll, int y_ll, int wdth, double t) {
   int x, y, count;
   double v, max_v;

   // Search the neighborhood for the maximum value
   count = 0;
   max_v = -1.0e30;
   for(x=x_ll; x<x_ll+wdth; x++) {
      if(x < 0 || x >= wd->get_nx()) continue;
   
      for(y=y_ll; y<y_ll+wdth; y++) {
         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         v = wd->get_xy_double(x, y);
         if(v > max_v) max_v = v;
         count++;
      } // end for y
   } // end for x
   
   // Check whether enough valid grid points were found to trust
   // the maximum value computed
   if( (double) count/(wdth*wdth) < t
    || count == 0) {
      max_v = bad_data_double;
   }
   
   return(max_v);
}

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

double interp_median(WrfData *wd, int x_ll, int y_ll, int wdth, double t) {
   double *data;
   int x, y, count;
   double v, median_v;
   
   // Allocate space to store the data points
   data = new double [wdth*wdth];
   
   // Search the neighborhood for valid data points
   count = 0;
   for(x=x_ll; x<x_ll+wdth; x++) {
      if(x < 0 || x >= wd->get_nx()) continue;
   
      for(y=y_ll; y<y_ll+wdth; y++) {
         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         v = wd->get_xy_double(x, y);
         data[count] = v;
         count++;
      } // end for y
   } // end for x
   
   // Check whether enough valid grid points were found to compute
   // a median value
   if( (double) count/(wdth*wdth) < t
    || count == 0) {
      median_v = bad_data_double;
   }
   else {
      sort(data, count);
      median_v = percentile(data, count, 0.50);
   }
   
   if(data) { delete [] data; data = (double *) 0; }
   
   return(median_v);
}

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

double interp_uw_mean(WrfData *wd, int x_ll, int y_ll, int wdth, double t) {
   int x, y, count;
   double v, sum, uw_mean_v;

   // Sum up the valid data in the neighborhood
   count = 0;
   sum = 0.0;
   for(x=x_ll; x<x_ll+wdth; x++) {
      if(x < 0 || x >= wd->get_nx()) continue;
   
      for(y=y_ll; y<y_ll+wdth; y++) {
         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         v = wd->get_xy_double(x, y);
         sum += v;
         count++;
      } // end for y
   } // end for x
   
   // Check whether enough valid grid points were found to compute 
   // a mean value
   if( (double) count/(wdth*wdth) < t
    || count == 0) {
      uw_mean_v = bad_data_double;
   }
   else {
      uw_mean_v = sum/count;
   }
   
   return(uw_mean_v);
}

////////////////////////////////////////////////////////////////////////
//
// Compute the distance-weighted mean using Shepards Method
//
////////////////////////////////////////////////////////////////////////

double interp_dw_mean(WrfData *wd, int x_ll, int y_ll, int wdth,
                      double obs_x, double obs_y, int i_pow, double t) {
   double *data, *dist;
   int i, x, y, count;
   double v, dw_mean_v, w, wght_sum;
   
   // Allocate space to store the data points
   data = new double [wdth*wdth];
   dist = new double [wdth*wdth];
   
   // Search the neighborhood for valid data points
   count = 0;
   for(x=x_ll; x<x_ll+wdth; x++) {
      if(x < 0 || x >= wd->get_nx()) continue;
   
      for(y=y_ll; y<y_ll+wdth; y++) {
         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         v = wd->get_xy_double(x, y);
         data[count] = v;
         dist[count] = sqrt(pow((obs_x-x), 2) + pow((obs_y-y), 2));
         count++;
      } // end for y
   } // end for x
   
   // Check whether enough valid grid points were found to compute
   // a distance weighted mean value
   if( (double) count/(wdth*wdth) < t
    || count == 0) {
      dw_mean_v = bad_data_double;
   }
   else {
   
      // Compute the sum of the weights
      for(i=0, wght_sum=0; i<count; i++) {
      
         // If the distance is very close to zero,
         // break out of the loop
         if(dist[i] <= 0.001) break;

         // Otherwise, compute the sum of the weights
         wght_sum += pow(dist[i], -1*i_pow);
      }
      
      for(i=0, dw_mean_v=0; i<count; i++) {
         
         // If the distance is very close to zero, set the dw_mean_v
         // value to the current data value and break out of the loop
         if(dist[i] <= 0.001) {
            dw_mean_v = data[i];
            break;
         }
         
         w = pow(dist[i], -1*i_pow)/wght_sum;
         dw_mean_v += w*data[i];
      }
   }
   
   if(data) { delete [] data; data = (double *) 0; }
   if(dist) { delete [] dist; dist = (double *) 0; }
   
   return(dw_mean_v);
}

////////////////////////////////////////////////////////////////////////
//
// Compute a least-squares fit interpolation
//
////////////////////////////////////////////////////////////////////////

double interp_ls_fit(WrfData *wd, int x_ll, int y_ll, int wdth, 
                     double obs_x, double obs_y, double t) {
   int i, j, x, y, count;
   const int N  = wdth;
   const int N2 = N*N;
   const double alpha     = (N2*(N2 - 1.0))/12.0;
   const double beta      = 0.5*(N - 1.0);
   const double x_center  = x_ll + beta;
   const double y_center  = y_ll + beta;
   const double u_test    = obs_x - x_center;
   const double v_test    = obs_y - y_center;
   double A, B, C;
   double suz, svz, sz;
   double u, v, z;

   if(N < 2) {
      cerr << "\n\nERROR: interp_ls_fit() -> "
           << "the interpolation width (" << N 
           << ") must be set >= 2\n\n" << flush;

      exit (1);
   }

   suz = svz = sz = 0.0;

   // Search the neighborhood
   count = 0;
   for(i=0; i<N; i++) {

      u = i - beta;
      x = x_ll + i;
      
      if(x < 0 || x >= wd->get_nx()) continue;

      for(j=0; j<N; j++) {

         v = j - beta;
         y = y_ll + j;

         if(y < 0 || y >= wd->get_ny()) continue;
         if(wd->is_bad_data(x, y))      continue;

         z = wd->get_xy_double(x, y);
         count++;

         suz += u*z;
         svz += v*z;
         sz  += z;
      }
   }

   A = suz/alpha;
   B = svz/alpha;
   C = sz/N2;

   z = A*u_test + B*v_test + C;
   
   // Check for not enough valid data
   if( (double) count/N2 < t || count == 0) {
      z = bad_data_double;
   }

   return(z);
}

///////////////////////////////////////////////////////////////////////////////
//
// Convert a field of data into its corresponding field of ranks of that data.
// Return the number of valid data points that were ranked and keep track
// of the number of rank ties.
//
///////////////////////////////////////////////////////////////////////////////

int compute_rank(const WrfData &wd, WrfData &wd_rank, double *data_rank, int &ties) {
   int x, y, n, i;
   double *data;
   int *data_loc;

   // Arrays to store the raw data values to be ranked, their locations, 
   // and their computed ranks.  The ranks are stored as doubles since
   // they can be set to 0.5 in the case of ties.
   data      = new double [wd.get_nx()*wd.get_ny()];
   data_loc  = new int    [wd.get_nx()*wd.get_ny()];
   
   // Search the input field for valid data and keep track of its location 
   n = 0;
   for(x=0; x<wd.get_nx(); x++) {
      for(y=0; y<wd.get_ny(); y++) {

         if(!wd.is_bad_data(x, y)) { 
            data[n] = wd.get_xy_double(x, y);
            data_loc[n] = wd.two_to_one(x, y);
            n++;
         }
      }
   }
   
   // Compute the rank of the data and store the ranks in the data_rank array
   // Keep track of the number of ties in the ranks
   ties = rank(data, data_rank, n);
   
   // Set up the wd_rank object
   wd_rank.set_size(wd.get_nx(), wd.get_ny());
   wd_rank.set_m((double) n/wrfdata_int_data_max);
   wd_rank.set_b(0.0);
   
   // Assign the ranks to the wd_rank field
   for(i=0; i<n; i++) {
      wd_rank.one_to_two(data_loc[i], x, y);
      wd_rank.put_xy_double(data_rank[i], x, y);
   }

   // Deallocate memory
   if(data)      { delete [] data;      data = (double *) 0; }
   if(data_loc)  { delete [] data_loc;  data_loc = (int *) 0; }

   return(n);   
}

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