[gradsusr] Plotting Hwind Analysis data in Grads

Rick Danielson Rick.Danielson at phys.ocean.dal.ca
Fri Aug 17 01:37:31 EDT 2012


Hi Sam,

      I'm not sure about GrADs reading ASCII directly (you might
need to write the grid as a binary file first, but perhaps another
GrADs user can help there).  My solution has a different challenge
if you're up to first installing NetCDF and udunits packages from
Unidata.  If so, the attached C program can then be compiled using
a command like this (all one line):

cc -o hwnd2nc hwnd2nc.c -L/contrib/netcdf-3.6.2/gcc4/lib -lnetcdf
-L/glade/home/rickd/soft/udunits/lib -ludunits2 -lm

(where my libnetcdf.a is installed in /contrib/netcdf-3.6.2/gcc4/lib
and libudunits2.a is in /glade/home/rickd/soft/udunits/lib and header
files in the attached hwnd2nc.c program would also need attention).
Executing "hwnd2nc AL122005_0828_1800 should yield an H*Wind data file
in NetCDF, which GrADs can then sdfopen (see other file attached).

Cheers,
Rick


> Hi all,
> 
> I'm wondering if anyone has any experience plotting Hwind analysis data in grads and, if so, could assist me in plotting it?
> 
> Attached is an example gridded data file for one time for hurricane Katrina.  The attached file contains useful information to create a ctl file for the
> data, such as the lat/lon bounds of the grid, but I'm a bit lost otherwise.
> 
> Thanks for any help you can provide.
> 
> Best Regards,
> Sam
-------------- next part --------------
/*
 * This program is designed to create a netcdf template (i.e. define the
 * fundamental dimensions, variables, and attributes, including data for
 * the coordinate variables) for surface-level spherical coordinate data.
 * The input dates are assumed to be in increasing order - RD November 1999.
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>
#include "/contrib/netcdf-3.6.2/gcc4/include/netcdf.h"
#include "/glade/home/rickd/soft/udunits/include/udunits.h"

#define LEN        999
#define LOTS       4000                      /* number of domain lats/lons (max) */
#define LEVELS     1                         /* number of pressure levels */
#define MISSING    -9e9                      /* generic missing value */

#define DIST       0                         /* distance from gridbox to cyclone center (km) */
#define LATS       1                         /* gridbox  latitude (deg) */
#define LONS       2                         /* gridbox longitude (deg) */
#define UWND       3                         /* zonal      wind component (m/s) */
#define VWND       4                         /* meridional wind component (m/s) */
#define PARAMS     5                         /* number of parameters */

void sarwrite(char file[], char param[], char date[], float *array, int arraylen);
void timematch(int ncid, char date[], int *match);
int get_mem2Dfloat(float ***array, int dimb, int dima);
int get_mem3Dfloat(float ****array, int dimc, int dimb, int dima);
void no_mem_exit(char *message);

float ***data, *array;

main(int argc, char *argv[])
{
    time_t time_val;                                                          /* time declarations (first!) */
    struct tm tma;
    char *month[12] = {"January","February","March","April","May","June","July",
                       "August","September","October","November","December"};

    FILE *fpa, *fopen();                                                      /* basic declarations */
    char fila[LEN], filb[LEN], filc[LEN];
    char date[LEN], datempa[LEN], datempb[LEN], line[5*LOTS];
    int a, b, c, d, count;

    int status, ncid;                                                         /* netCDF declarations */
    int time_dim, level_dim, lat_dim, lon_dim;
    int time_id,  level_id,  lat_id,  lon_id,  par_id[PARAMS], dims[4];
    short short_val;
    float float_val, float_range[2];
    double double_range[2];
    char attrib[LEN], lina[LEN], linb[LEN];
    static size_t begin[3] = {0,0,0}, level_end = LEVELS, lat_end, lon_end;
    size_t time_end, indic_end[3];

    char udunita[LEN], udunitb[LEN];                                          /* udunits declarations */
    utUnit utunita, utunitb;
    double slope, intercept, times[LOTS];
    char predate[] = "hours since ";
    char postdate[] = ":00:0.0";

    char *long_name[PARAMS] = {"Distance from cyclone center",
                               "Latitude of wind vector cell",
                               "Longitude of wind vector cell",
                               "Zonal wind component",
                               "Meridional wind component"};
    char *parname[PARAMS] = {    "dist", "lats", "lons", "uwnd", "vwnd"};
    char *units[PARAMS] = {        "km",  "deg",  "deg",  "m/s",  "m/s"};
    float add_offset[PARAMS] = {     0.,     0.,     0.,     0.,     0.};
    float scale_factor[PARAMS] = {   1.,     1.,     1.,     1.,     1.};
    float missing_value[PARAMS] = {-9e9,   -9e9,   -9e9,   -9e9,   -9e9};
    int numdims[PARAMS] = {           3,      3,      3,      3,      3};
    float lat[LOTS], lon[LOTS], dis[LOTS];

    if (argc != 2) {
      printf("Usage: %s input_file\n",argv[0]);
      printf(" e.g.: %s AL132003_0913_2204marine\n",argv[0]);
      exit(1);
    }
    strcpy(fila,argv[1]);                                                     /* get the date from the input file name */
    date[ 0] = fila[ 4];
    date[ 1] = fila[ 5];
    date[ 2] = fila[ 6];
    date[ 3] = fila[ 7];
    date[ 4] = '-';
    date[ 5] = fila[ 9];
    date[ 6] = fila[10];
    date[ 7] = '-';
    date[ 8] = fila[11];
    date[ 9] = fila[12];
    date[10] = '-';
    date[11] = fila[14];
    date[12] = fila[15];
    date[13] = '\0';

    if ((fpa = fopen(fila,"r")) == NULL) {                                    /* read grid dimensions from the input file */
      fprintf(stderr,"ERROR : couldn't open %s\n",fila);                      /* (with lat order reversed: from north to south) */
      exit(1);
    }
    printf("reading %s\n",fila);

    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa);
    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa) ; fscanf(fpa,"%d",&lon_end);
    if (lon_end > LOTS) {printf("ERROR : lon_end %d > LOTS %d\n",lon_end,LOTS) ; exit(1);}
    for (a = 0; a < lon_end; a++) fscanf(fpa,"%f",&lon[a]);
/* for (a = 0; a < lon_end; a++) printf("%f ",lon[a]); */
    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa) ; fscanf(fpa,"%d",&lat_end);
    if (lat_end > LOTS) {printf("ERROR : lat_end %d > LOTS %d\n",lat_end,LOTS) ; exit(1);}
    for (a = 0; a < lat_end; a++) fscanf(fpa,"%f",&lat[lat_end-a-1]);

    get_mem3Dfloat(&data,PARAMS,lat_end,lon_end);                             /* allocate space for grids and calculate the */
    array = (float *)calloc(lat_end*lon_end,sizeof(float));                   /* gridded distance to the cyclone center */
    for (a = 0; a < lat_end; a++)
      for (b = 0; b < lon_end; b++)
        data[DIST][a][b] = powf(lat[a]*lat[a]+lon[b]*lon[b],0.5);

    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa) ; fscanf(fpa,"%d",&lon_end);
    if (lon_end > LOTS) {printf("ERROR : lon_end %d > LOTS %d\n",lon_end,LOTS) ; exit(1);}
    for (a = 0; a < lon_end; a++) fscanf(fpa,"%f",&lon[a]);
    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa) ; fscanf(fpa,"%d",&lat_end);
    if (lat_end > LOTS) {printf("ERROR : lat_end %d > LOTS %d\n",lat_end,LOTS) ; exit(1);}
    for (a = 0; a < lat_end; a++) fscanf(fpa,"%f",&lat[lat_end-a-1]);

    for (a = 0; a < lat_end; a++)                                             /* and get lats and lons (with lat order reversed) */
      for (b = 0; b < lon_end; b++) {
        data[LATS][a][b] = lat[a];
        data[LONS][a][b] = lon[b];
      }

    fgets(line,LEN,fpa) ; fgets(line,LEN,fpa) ; fscanf(fpa,"%*s %*s");
/* printf("line:%s",line);
printf("lon_end:%d\n",lon_end);
printf("lat_end:%d\n",lat_end); */
    for (a = 0; a < lat_end; a++)
      for (b = 0; b < lon_end; b++)
        fscanf(fpa,"%*s %f, %f)",&data[UWND][lat_end-a-1][b],&data[VWND][lat_end-a-1][b]);

    a = 0;
    if (utInit("") != 0) {
      printf("couldn't initialize Unidata units library\n");
      exit(2);
    }
    strcpy(udunita,predate);                                                  /* form the desired date */
    strcat(udunita,date);
    strcat(udunita,postdate);
    udunita[22] = ' ';
    strcpy(udunitb,"hours since 1-1-1 00:00:0.0");                            /* and the reference date */
    if (utScan(udunita,&utunita) != 0 || utScan(udunitb,&utunitb) != 0)       /* to obtain the number of hours */
      printf("either %s or %s is incomprehensible\n",udunita,udunitb);        /* from the reference date */
    if (utConvert(&utunita,&utunitb,&slope,&intercept) != 0)                  /* to the desired date */
      printf("couldn't convert %s to %s\n",udunita,udunitb);
    times[a++] = intercept;
    utTerm();
    count = a;

    strcpy(filb,argv[1]) ; strcat(filb,".nc");                                /* form the template file name */
    status = nc_create(filb, NC_CLOBBER, &ncid);                              /* create file and enter define mode */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    printf("writing %s\n",filb);

    status = nc_def_dim(ncid, "time", NC_UNLIMITED, &time_dim);               /* define dimensions */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_def_dim(ncid, "lat", lat_end, &lat_dim);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_def_dim(ncid, "lon", lon_end, &lon_dim);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    printf("with dimensions of time %d, lat %d, and lon %d\n",count,lat_end,lon_end);

    dims[0] = time_dim;                                                       /* define dimension variables */
    status = nc_def_var(ncid, "time", NC_DOUBLE, 1, dims, &time_id);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    dims[0] = lat_dim;
    status = nc_def_var(ncid, "lat", NC_FLOAT, 1, dims, &lat_id);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    dims[0] = lon_dim;
    status = nc_def_var(ncid, "lon", NC_FLOAT, 1, dims, &lon_id);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    printf("including variables");
    dims[0] = time_dim;                                                       /* define variables */
    for (a = 0; a < PARAMS; a++) {
      if (numdims[a] == 3) {
        dims[1] = lat_dim;
        dims[2] = lon_dim;
      }
      else {
        dims[1] = level_dim;
        dims[2] = lat_dim;
        dims[3] = lon_dim;
      }
      status = nc_def_var(ncid, parname[a], NC_FLOAT, numdims[a], dims, &par_id[a]);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      printf(" %s",parname[a]);
    }
    printf("\n");

    strcpy(attrib,"long_name");                                               /* assign time attributes */
    strcpy(line,"Time");
    status = nc_put_att_text(ncid, time_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"units");
    strcpy(line,"hours since 1-1-1 00:00:0.0");
    status = nc_put_att_text(ncid, time_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"actual_range");
    double_range[0] = times[0];
    double_range[1] = times[count-1];
    status = nc_put_att_double(ncid, time_id, attrib, NC_DOUBLE, 2, double_range);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"delta_t");
    strcpy(line,"0000-00-00 06:00:00");
    status = nc_put_att_text(ncid, time_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    strcpy(attrib,"long_name");                                               /* assign lat attributes */
    strcpy(line,"Latitude");
    status = nc_put_att_text(ncid, lat_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"units");
    strcpy(line,"degrees_north");
    status = nc_put_att_text(ncid, lat_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"actual_range");
    float_range[0] = lat[0];
    float_range[1] = lat[lat_end-1];
    status = nc_put_att_float(ncid, lat_id, attrib, NC_FLOAT, 2, float_range);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    strcpy(attrib,"long_name");                                               /* assign lon attributes */
    strcpy(line,"Longitude");
    status = nc_put_att_text(ncid, lon_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"units");
    strcpy(line,"degrees_east");
    status = nc_put_att_text(ncid, lon_id, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"actual_range");
    float_range[0] = lon[0];
    float_range[1] = lon[lon_end-1];
    status = nc_put_att_float(ncid, lon_id, attrib, NC_FLOAT, 2, float_range);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    for (a = 0; a < PARAMS; a++) {                                            /* assign parameter attributes */
      strcpy(attrib,"long_name");
      strcpy(line,long_name[a]);
      status = nc_put_att_text(ncid, par_id[a], attrib, strlen(line)+1, line);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      strcpy(attrib,"units");
      strcpy(line,units[a]);
      status = nc_put_att_text(ncid, par_id[a], attrib, strlen(line)+1, line);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      strcpy(attrib,"add_offset");
      float_val = add_offset[a];
      status = nc_put_att_float(ncid, par_id[a], attrib, NC_FLOAT, 1, &float_val);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      strcpy(attrib,"scale_factor");
      float_val = scale_factor[a];
      status = nc_put_att_float(ncid, par_id[a], attrib, NC_FLOAT, 1, &float_val);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      strcpy(attrib,"missing_value");
      float_val = missing_value[a];
      status = nc_put_att_float(ncid, par_id[a], attrib, NC_FLOAT, 1, &float_val);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    }

    strcpy(attrib,"Conventions");                                             /* assign global attributes */
    strcpy(line,"COARDS");
    status = nc_put_att_text(ncid, NC_GLOBAL, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"title");
    sprintf(line,"H*Wind data for %s",argv[1]);
    status = nc_put_att_text(ncid, NC_GLOBAL, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"history");
    time_val = time(NULL);
    tma = *localtime(&time_val);
    strcpy(line,"H*Wind analyses can be obtained from storm.aoml.noaa.gov/\n");
    sprintf(lina,"hwind/Analysis/ - RD %s %d.",month[tma.tm_mon],1900+tma.tm_year);
    strcat(line,lina);
    status = nc_put_att_text(ncid, NC_GLOBAL, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    strcpy(attrib,"description");
    strcpy(line,"Analyses were performed using conventional settings with\n");
    strcat(line,"observations inside a 6h centered window (usually) considered.");
    status = nc_put_att_text(ncid, NC_GLOBAL, attrib, strlen(line)+1, line);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    status = nc_enddef(ncid);                                                 /* leave define mode */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    time_end = (long)count;                                                   /* store dimension variables */
    status = nc_put_vara_double(ncid, time_id, begin, &time_end, times);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_put_vara_float(ncid, lat_id, begin, &lat_end, lat);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_put_vara_float(ncid, lon_id, begin, &lon_end, lon);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    status = nc_close(ncid);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    d = 0 ; for (a = 0; a < lat_end; a++) for (b = 0; b < lon_end; b++) array[d++] = data[DIST][a][b];
    sarwrite(filb,"dist",date,array,lat_end*lon_end);
    d = 0 ; for (a = 0; a < lat_end; a++) for (b = 0; b < lon_end; b++) array[d++] = data[LATS][a][b];
    sarwrite(filb,"lats",date,array,lat_end*lon_end);
    d = 0 ; for (a = 0; a < lat_end; a++) for (b = 0; b < lon_end; b++) array[d++] = data[LONS][a][b];
    sarwrite(filb,"lons",date,array,lat_end*lon_end);
    d = 0 ; for (a = 0; a < lat_end; a++) for (b = 0; b < lon_end; b++) array[d++] = data[UWND][a][b];
    sarwrite(filb,"uwnd",date,array,lat_end*lon_end);
    d = 0 ; for (a = 0; a < lat_end; a++) for (b = 0; b < lon_end; b++) array[d++] = data[VWND][a][b];
    sarwrite(filb,"vwnd",date,array,lat_end*lon_end);
    exit(0);
}

void sarwrite(char file[], char param[], char date[], float *array, int arraylen)
{
    int match, loopa, loopb, datalen, mesgflag;
    double packedval;
    int status, ncid, dim_id, par_id, par_ndims, par_dims[NC_MAX_VAR_DIMS], par_natts;
    char dim_name[NC_MAX_NAME], dim_specs[NC_MAX_VAR_DIMS*NC_MAX_NAME];
    nc_type par_type;
    float offset, scale, badval;
    size_t dim_len, origin[NC_MAX_VAR_DIMS], interval[NC_MAX_VAR_DIMS];

    if (strcmp(param,"") == 0)                                                /* just return if the */
      return;                                                                 /* parameter name is null */

    status = nc_open(file, NC_WRITE, &ncid);                                  /* open the data file */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    status = nc_inq_varid(ncid, param, &par_id);                              /* get parameter information */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_inq_var(ncid, par_id, 0, &par_type, &par_ndims, par_dims, &par_natts);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    datalen = 1;                                                              /* get dimension information */
    strcpy(dim_specs,"");
    for (loopa = 0; loopa < par_ndims; loopa++) {
      status = nc_inq_dim(ncid, par_dims[loopa], dim_name, &dim_len);
      if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
      strcat(dim_specs,dim_name);
      strcat(dim_specs,"=");
      if (strcmp(dim_name,"time") == 0) {
        timematch(ncid, date, &match);
        origin[loopa] = match;
        interval[loopa] = 1;
        sprintf(dim_name,"%d",match+1);
        strcat(dim_specs,dim_name);
      }
      else {
        origin[loopa] = 0;
        interval[loopa] = dim_len;
        datalen *= (int)dim_len;
        sprintf(dim_name,"1-%d",(int)dim_len);
        strcat(dim_specs,dim_name);
      }
      if (loopa != par_ndims-1)
        strcat(dim_specs," ");
    }

    if (arraylen < datalen) {                                                 /* exit if the number of data values */
      printf("ERROR: ");                                                      /* to write is less than required */
      printf("%d %s values available to fill in %s ",datalen,param,file);
      printf("but array has only %d values\n",arraylen);
      exit(-1);
    }
    else if (arraylen > datalen) {                                            /* warn if the number of data values */
      printf("WARNING: ");                                                    /* to write is more than required */
      printf("%d %s values available to fill in %s ",datalen,param,file);
      printf("but array has %d values\n",arraylen);
    }

    mesgflag = 0;
    for (loopa = 0; loopa < datalen; loopa++) {                               /* check the grid */
      if (array[loopa] == MISSING) {                                          /* and if any data values are */
/*      array[loopa] = MISSING;                                                * missing then reset these according */
        mesgflag++;                                                           /* to the internal convention (for short) */
/*      badval = array[loopa];
      }
      else
        array[loopa] = array[loopa];  */}
    }

    if (mesgflag == 1) {                                                      /* and notify the user */
      printf("WARNING: ");
      printf("one %s value (%f) exceeds the valid packing range ",param,badval);
      printf("so the corresponding packed value has been reset to MISSING (%f)\n",MISSING);
    }
    else if (mesgflag > 1) {
      printf("WARNING: ");
      printf("%d %s values (e.g. %f) exceed the valid packing range ",mesgflag,param,badval);
      printf("so all corresponding packed values have been reset to MISSING (%f)\n",MISSING);
    }

    printf("writing %s %s %-12s %s\n",file,date,param,dim_specs);             /* finally write the packed grid */
    status = nc_put_vara_float(ncid, par_id, origin, interval, array);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}

    status = nc_close(ncid);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
}

#define MATCHTIMES 10000
#define MATCHLEN   500

void timematch(int ncid, char date[], int *match)
{
    size_t count, timlen;                                                     /* netcdf declarations */
    int status, par_id, dim_id;
    double times[MATCHTIMES];

    char unita[MATCHLEN], unitb[MATCHLEN];                                    /* udunits declarations */
    utUnit utunita, utunitb;
    double slope, intercept;

    status = nc_inq_dimid(ncid, "time", &dim_id);                             /* count and retrieve all */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}    /* (referenced) date values */
    status = nc_inq_dimlen(ncid, dim_id, &count);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_inq_varid(ncid, "time", &par_id);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_get_var_double(ncid, par_id, times);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_inq_attlen(ncid, par_id, "units", &timlen);
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    status = nc_get_att_text(ncid, par_id, "units", unitb);                   /* get the reference date */
    if (status != NC_NOERR) {printf("%s\n",nc_strerror(status));exit(-1);}
    unitb[timlen] = '\0';

    if (utInit("") != 0) {
      printf("couldn't initialize Unidata units library\n");
      exit(2);
    }

    sprintf(unita,"hours since %s:00:0.0", date);                             /* find the referenced date value */
    unita[22] = ' ';                                                          /* for the desired date */
    utScan(unita,&utunita);
    utScan(unitb,&utunitb);
    utConvert(&utunita,&utunitb,&slope,&intercept);
    utTerm();

    for (*match = 0; *match < count; (*match)++)                              /* match the desired date value with */
      if (intercept == times[*match])                                         /* one from the input netcdf file */
        break;

    if (*match == count) {                                                    /* but if no match was found then */
      printf("ERROR: ");                                                      /* print a warning and exit */
      printf("could not match %s with any date in the netCDF file\n",date);
      exit(-1);
    }
}

int get_mem2Dfloat(float ***array, int dimb, int dima)
{
    int i;

    if (( *array     = (float**)calloc(dimb,     sizeof(float*))) == NULL)
      no_mem_exit("ERROR allocating memory in get_mem2Dfloat");
    if (((*array)[0] = (float* )calloc(dimb*dima,sizeof(float ))) == NULL)
      no_mem_exit("ERROR allocating memory in get_mem2Dfloat");
    for (i = 1; i < dimb; i++)
      (*array)[i] =  (*array)[i-1] + dima;

    return dimb*dima*sizeof(float);
}

int get_mem3Dfloat(float ****array, int dimc, int dimb, int dima)
{
    int  i;

    if (((*array) = (float***)calloc(dimc,sizeof(float**))) == NULL)
      no_mem_exit("ERROR allocating memory in get_mem3Dfloat");
    for (i = 0; i < dimc; i++)
      get_mem2Dfloat((*array)+i, dimb, dima);

    return dimc*dimb*dima*sizeof(float);
}

void no_mem_exit(char *message)
{
    fprintf(stderr,"\n%s\n\n",message);
    exit(1);
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: AL122005_0828_1800.nc
Type: application/x-netcdf
Size: 521652 bytes
Desc: 
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20120817/80265024/attachment-0003.nc 


More information about the gradsusr mailing list