/* reads a MAPS netCDF file and generates sounding data
 * 30-Sep-98 - outputs CAPE, CIn, and Helic
 *           - changed to Marks NEW FSL sounding format
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "netcdf.h"
#define MAX_LEVELS 200


main (int argc, char *argv[])
{
  clock_t tstart, tend;
  double elapsed;
  int n_reads=0;

  long start[]={0,0,0,0,0};
  long count[]={1,1,1,1,1};
  long stride[]={1,1,1,1,1};
  float VPT[MAX_LEVELS];
  char *filename;
  char *source;  /* MAPS or RUC */
  char *fieldname;
  nc_type data_type;
  int i,data_id,nc_id,data_ndims,tot_len,startx,starty;
  int data_dimids[5];
  long data_dim_values[5];
  int icount;			/* counts output levels */
  int p_levels;
 
  source = "RUC";
  filename = "/w3/rapb/ruc_native_40.nc";
  n_reads = atoi(argv[1]);

  tstart = clock();

  for(icount=0;icount<n_reads;icount++) {
  ncopts=NC_VERBOSE;			/* don't die on read errors */
  nc_id=ncopen(filename,NC_NOWRITE);
  if(nc_id == -1) {
    fprintf(stderr,"!!No data found for %s\n",filename);
    exit(2);
  }

  /* now read the  file. */
  fieldname="u";
  data_id = ncvarid(nc_id,fieldname);
  
  ncvarinq(nc_id,data_id,(char *) 0, &data_type, &data_ndims,
	   data_dimids, (int *) 0);
 
  /*
  fprintf(stderr,"Variable %s is of type %d, with %d dimensions.\n"
	 ,fieldname,data_type,data_ndims,
	 tot_len);
  */

  /* find the values of each dimension */
  for (i=0;i<data_ndims;i++) {
    ncdiminq(nc_id,data_dimids[i],(char *) 0,&data_dim_values[i]);
    /* fprintf(stderr,"%d: dimension %d has value %d\n",i,data_dimids[i],
	   data_dim_values[i]);*/
  }

  /* get x,y point */
  startx = 50;
  starty = 40;
  start[2]=starty;
  start[3]=startx;

  p_levels = data_dim_values[1];
  count[1]=p_levels;
  stride[2]=data_dim_values[2];
  stride[3]=data_dim_values[3];
/* fprintf(stderr,"%d %d %d %d\n",stride[0],stride[1],stride[2],stride[3]); */

 /* read one value from each level using a generalized hyperslab   */
  data_id=ncvarid(nc_id,"vpt");
  ncvargetg(nc_id,data_id,start,count,stride,(char *)0,(void *)VPT);

  /*
  icount = 0;		
   for(i=0;i<p_levels;i++) {
     printf("%d: t %f\n",i,VPT[i]);
   }
  */
  ncclose(nc_id);
  }
  tend = clock();
  elapsed = ((double) (tend - tstart)) / CLOCKS_PER_SEC;
  printf("C: elapsed time for %d reads is %f seconds\n",n_reads,elapsed);
}
