mm5-nc2F5.c

Converts MM5 files, given in NetCDF format, to F5, as time-dependent 3D unidistant scalar and vector fields.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <netcdf.h>
#include <hdf5.h>
#include <F5/F5uniform.h>
#include <F5/F5coordinates.h>
#include <F5/F5X.h>
 
 
#ifdef DEBUG_ON
# define DEBUG(fmt, ...) printf("DEBUG: " fmt "\n", __VA_ARGS__)
#else
# define DEBUG(...)
#endif
 
#define MAXDIMS (NC_MAX_VAR_DIMS)
#define MAXGROUPS 10
#define MAXPATH 1024
#define MAXATTR 256
#define VARSEP ",; "
 
 
struct varinfo
{
  char *name;
  nc_type type;
  int ndims, id;
  size_t *dim;
  int *ivar, size;
};
 
struct vardata
{
  struct varinfo *ncref;
  union
  {
    void *vptr;
    double *dptr;
    float *fptr;
    int *iptr;
    short *sptr;
  };
  hid_t h5type;
};
 
 
typedef struct varinfo vinfo;
typedef struct vardata vdata;
 
 
void error(char *msg)
{
  fprintf(stderr, "ERROR: %s\n", msg);
  exit(1);
}
 
vinfo *get_varinfo(int fd, int id)
{
  int i, dimid[MAXDIMS], natt;
  char nm[NC_MAX_NAME+1];
  vinfo *v = 0;
  
  v = malloc(sizeof(vinfo));
  if (nc_inq_var(fd, id, nm, &(v->type), &(v->ndims), dimid,
                 &natt) != NC_NOERR)
    error("cannot retrieve variable information");
 
  v->id = id;
  strcpy(v->name = malloc(strlen(nm)+1), nm);
  if (v->ndims > 0)
  {
    int tot = 1;
 
    v->ivar = malloc(v->ndims*sizeof(int));
    v->dim = malloc(v->ndims*sizeof(size_t));
    for (i = 0; i < v->ndims; i++)
    {
      v->ivar[i] = dimid[i];
      if (nc_inq_dimlen(fd, dimid[i], v->dim+i) != NC_NOERR)
        error("failed to retrieve dimension length");
      DEBUG("dimlen of index id=%d is %d", dimid[i], v->dim[i]);
      tot *= v->dim[i];
    }
    v->size = tot;
  }
  else
  {
    v->ivar = v->dim = 0;
    v->size = 1;
  }
  DEBUG("variable '%s' info: id=%d ndims=%d size=%d", nm, id, v->ndims, v->size);
 
  return v;
}
 
vdata **get_ncvars(int fd, int nvars, char *vars[])
{
  int i, vid, nd;
  size_t start[MAXDIMS];
  nc_type type;
  vdata **outvar = malloc(sizeof(vdata)*nvars);
 
  memset(start, 0, sizeof(size_t)*MAXDIMS);
 
  for (i = 0; i < nvars; i++)
  {
    int rc, k;
    vinfo *vi;
    vdata *vd = malloc(sizeof(vdata));
 
    outvar[i] = vd;
    if (nc_inq_varid(fd, vars[i], &vid) != NC_NOERR)
      error("cannot find requested variable");
    vd->ncref = vi = get_varinfo(fd, vid);
  }
 
  return outvar;
}
 
void read_ncvars(int fd, vinfo *ui, vdata **vd, int nvars)
{
  int i;
  size_t start[MAXDIMS];
 
  memset(start, 0, sizeof(size_t)*MAXDIMS);
 
  for (i = 0; i < nvars; i++)
  {
    vinfo *vi = vd[i]->ncref;
 
    switch (vi->type)
    {
    case NC_FLOAT:
      vd[i]->fptr = malloc(ui->size*sizeof(float));
      DEBUG("allocated %d bytes for \"%s\" data", ui->size*sizeof(float), vi->name);
      if (nc_get_vara_float(fd, vi->id, start, ui->dim, vd[i]->fptr) != NC_NOERR)
        error("failed to read float array");
      vd[i]->h5type = H5T_NATIVE_FLOAT;
      break;
    case NC_DOUBLE:
      vd[i]->dptr = malloc(ui->size*sizeof(double));
      DEBUG("allocated %d bytes for \"%s\" data", ui->size*sizeof(double), vi->name);
      if (nc_get_vara_double(fd, vi->id, start, ui->dim, vd[i]->dptr) != NC_NOERR)
        error("failed to read double array");
      vd[i]->h5type = H5T_NATIVE_DOUBLE;
      break;
    case NC_INT:
      vd[i]->iptr = malloc(ui->size*sizeof(int));
      DEBUG("allocated %d bytes for \"%s\" data", ui->size*sizeof(int), vi->name);
      if (nc_get_vara_int(fd, vi->id, start, ui->dim, vd[i]->iptr) != NC_NOERR)
        error("failed to read int array");
      vd[i]->h5type = H5T_NATIVE_INT;
      break;
    case NC_SHORT:
      vd[i]->sptr = malloc(ui->size*sizeof(short));
      DEBUG("allocated %d bytes for \"%s\" data", ui->size*sizeof(short), vi->name);
      if (nc_get_vara_short(fd, vi->id, start, ui->dim, vd[i]->sptr) != NC_NOERR)
        error("failed to read short array");
      vd[i]->h5type = H5T_NATIVE_SHORT;
      break;
    default:
      /* TODO: add other netcdf types */
      error("unrecognized variable type");
    }
  }
}
 
vinfo *get_time(int fd, char *tname)
{
  int tid;
  vinfo *tv;
 
  if (nc_inq_varid(fd, tname, &tid) != NC_NOERR)
    error("couldn't find time variable");
  if (!(tv = get_varinfo(fd, tid)))
    error("failed to extract time variable info");
  if ((tv->ndims > 1) || ((tv->ndims == 1) && (tv->dim[0] > 1)))
    error("multiple time steps within netCDF file are not supported yet");
  return tv;
}
 
vinfo *unify(vinfo *tinfo, vdata **vd, int vnum)
{
  int i, k;
  vinfo *vi, *m = malloc(sizeof(vinfo));
 
  for (i = 0; i < vnum; i++)
  {
    vi = vd[i]->ncref;
    if (!i)
    {
      m->dim = malloc(sizeof(size_t)*(m->ndims = vi->ndims));
      memcpy(m->dim, vi->dim, m->ndims*sizeof(size_t));
    }
    else
    {
      if (m->ndims != vi->ndims) return 0;
      for (k = 0; k < vi->ndims; k++)
      { /* simplistic for now: just chop of any additional grid points */
        if (m->dim[k] > vi->dim[k]) m->dim[k] = vi->dim[k];
      }
    }
  }
 
  for (m->size = 1, k = 0; k < m->ndims; k++)
  {
    DEBUG("unified dim[%d]=%d", k, m->dim[k]);
    m->size *= m->dim[k];
  }
 
  return m;
}
 
void dimopt(vinfo *ui)
{
  int k, l = 0;
 
  for (k = 0; k < ui->ndims; k++)
  {
    if (ui->dim[k] != 1) ui->dim[l++] = ui->dim[k];
    else
    {
      DEBUG("eliminated index %d", k);
    }
  }
  ui->ndims = l;
  DEBUG("number of unified dimensions set to %d:", l);
  for (k = 0; k < ui->ndims; k++)
    DEBUG("dim[%d]=%d", k, ui->dim[k]);
}
 
void usage(char *pname, int rc)
{
  fprintf(stderr, "Usage: %s [options] netcdf-files...\n", pname);
  fprintf(stderr, "options:\n"
          "  -v  var1[,var2]...   : variable group of interest (required)\n"
          "  -n  name             : dataset name for the preceding variable group\n"
          "  -t  name             : name of time variable\n"
          "  -o  pathname         : output file name (default: out.f5)\n"
          "  -r  float            : domain resolution (meters)\n");
  exit(rc);
}
 
char **optvars(char *s, int *n)
{
  int space = MAXDIMS, done = 0;
  char *p, **tab = malloc(sizeof(char *)*space);
 
  for (*n = 0; s[0]; (*n)++)
  {
    if (*n >= space)
    {
      space += MAXDIMS;
      tab = realloc(tab, sizeof(char *)*space);
    }
    p = strpbrk(s, VARSEP);
    if (!p) {p = s+strlen(s); done = 1;}
    tab[*n] = strncpy(malloc(p-s+1), s, p-s);
    tab[*n][p-s] = 0;
    s += p-s+(done == 0);
  }
 
  return tab;
}
 
F5Path *writef5_uni3d_scalar(hid_t f, double t, char *name,
                             F5_vec3_point_t *orig, F5_vec3_float_t *space,
                             vinfo *ui, vdata *vd)
{
  char tname[MAXATTR];
  hsize_t dim[3];
  F5Path *f5p;
 
  if (ui->ndims != 3) error("not a 3D grid");
 
  dim[2] = ui->dim[0];
  dim[1] = ui->dim[1];
  dim[0] = ui->dim[2];
 
  space->x /= ui->dim[0];
  space->y /= ui->dim[1];
  space->z /= ui->dim[2];
 
  DEBUG("scalar type %c= float",
        (H5Tequal(H5T_NATIVE_FLOAT, vd->h5type) > 0? '=': '!'));
  DEBUG("invoking F5Fwrite_uniform_cartesian3D() with dim=[%lld,%lld,%lld] name=\"%s\" data=%p", 
        dim[0], dim[1], dim[2], name, vd->vptr);
  f5p = F5Fwrite_uniform_cartesian3D(f, t, "UniformCartesian", orig, space,
                                     dim, name, vd->h5type, vd->vptr,
                                     0, H5P_DEFAULT);
  F5close(f5p);
 
  return f5p;
}
 
void writef5_uni3d_f3(hid_t f, double t, char *name,
                      F5_vec3_point_t *orig, F5_vec3_float_t *space,
                      vinfo *ui, vdata **vd)
{
  char tname[MAXATTR];
  hsize_t dim[3];
  const void *data[3];
  const float *fdata[3];
  F5Path *f5p;
  int i, rc, base;
  int NumOfDataValues;
  F5_vec3f_t* VectorData;
 
  if (ui->ndims != 3) error("not a 3D grid");
 
  dim[2] = ui->dim[0];
  dim[1] = ui->dim[1];
  dim[0] = ui->dim[2];
 
  data[0] = vd[0]->vptr;
  data[1] = vd[1]->vptr;
  data[2] = vd[2]->vptr;
 
  fdata[0] = (float*)vd[0]->vptr;
  fdata[1] = (float*)vd[1]->vptr;
  fdata[2] = (float*)vd[2]->vptr;
  /*Added by shalini to write array of vectors*/
  NumOfDataValues = dim[2]*dim[1]*dim[0];
  VectorData = malloc(sizeof(F5_vec3f_t)*NumOfDataValues);
  /*Convert each component array to the array of vectors*/
  for (i=0;i<NumOfDataValues;i++) {
        VectorData[i].x = fdata[0][i]; 
        VectorData[i].y = fdata[1][i]; 
        VectorData[i].z = fdata[2][i]; 
        /*fprintf(f,"Converted vector for base %d - {%f %f %f}\n",base, data[base],data[base+1],data[base+2]);*/
}
  if (!(f5p = F5Rcreate_uniform_cartesian3D(f, t, "UniformCartesian",
                                            orig, space, dim,
                                            0)))
    error("failed to create uniform grid");
 
  DEBUG("invoking F5Fwrites(%p, %s, 3, [%lld,%lld,%lld], float-3D, float-3D, "
        "[%p,%p,%p], default)", 
        f5p, name, dim[0], dim[1], dim[2], data[0], data[1], data[2]);
//  F5Fwrites(f5p, name, 3, dim, F5T_VEC3_FLOAT, F5T_VEC3_FLOAT, data,
        //    H5P_DEFAULT);
//Shalin- to write the tor  values as is ie array of structs
         F5Fwrite(f5p, name, /* fieldname */
                                3, dim, /* dataspace */
                                F5T_VEC3_FLOAT, /* type as to appear in the file */
                                F5T_VEC3_FLOAT, /* type as it resides in memory */
                                VectorData,     /* dataPtr */
                                F5P_DEFAULT);   /* property_id */
 
}
 
#if 0
void writef5_uni3d_f3(hid_t f, double t, char *name,
                      F5_vec3_point_t *orig, F5_vec3_float_t *space,
                      vinfo *ui, vdata **vd)
{
  char tname[MAXATTR];
  hsize_t dim[3];
  F5_vec3f_t *data;
  /*const void *data[3];*/
  F5Path *f5p;
  int i, rc;
 
  if (ui->ndims != 3) error("not a 3D grid");
 
  dim[0] = ui->dim[0];
  dim[1] = ui->dim[1];
  dim[2] = ui->dim[2];
 
  /*
  data[0] = vd[0]->vptr;
  data[1] = vd[1]->vptr;
  data[2] = vd[2]->vptr;
  */
 
  data = malloc(sizeof(float)*3*ui->size);
 
  for (i = 0; i < ui->size; i++)
  {
    data[i].x = vd[0]->fptr[i];
    data[i].y = vd[1]->fptr[i];
    data[i].z = vd[2]->fptr[i];
  }
 
  if (!(f5p = F5Rcreate_uniform_cartesian3D(f, t, "UniformCartesian",
                                            orig, space, dim,
                                            0)))
    error("failed to create uniform grid");
 
  DEBUG("invoking F5Fwrite(%p, %s, 3, [%lld,%lld,%lld], float-3D, float-3D, %p, default)", 
        f5p, name, dim[0], dim[1], dim[2], data);
  F5Fwrite(f5p, name, 3, dim, F5T_VEC3_FLOAT, F5T_VEC3_FLOAT, data,
           H5P_DEFAULT);
 
  /*
  DEBUG("invoking F5Fwrites(%p, %s, 3, [%lld,%lld,%lld], float-3D, float-3D, "
        "[%p,%p,%p], default)", 
        f5p, name, dim[0], dim[1], dim[2], data[0], data[1], data[2]);
  F5Fwrites(f5p, name, 3, dim, F5T_VEC3_FLOAT, F5T_VEC3_FLOAT, data,
            H5P_DEFAULT);
  */
}
 
void writef5_cart3d(hid_t f, double t, char *name,
                    F5_vec3_point_t *orig, F5_vec3_float_t *space,
                    vinfo *ui, vdata **vd, int nvars)
{
  char tpath[MAXATTR];
  hsize_t dim[3], tsz = 0;
  void *data[MAXDIMS];
  hid_t h5t;
  int i;
 
  dim[0] = ui->dim[0];
  dim[1] = ui->dim[1];
  dim[2] = ui->dim[2];
 
  snprintf(tpath, MAXATTR, "%s_type", name);
  h5t = H5Topen(f, tpath);
  if (h5t < 0)
  {
    tsz = H5Tget_size(vd[0]->h5type);
    h5t = H5Tcreate(H5T_COMPOUND, tsz*nvars);
    DEBUG("created compound type of %d bytes", tsz*nvars);
 
    for (i = 0; i < nvars; i++)
    {
      if (tsz != H5Tget_size(vd[i]->h5type))
        error("different sizes of vector component types");
      H5Tinsert(h5t, vd[i]->ncref->name, tsz*i, vd[i]->h5type);
      DEBUG("inserted \"%s\" member at offset %d", vd[i]->ncref->name, tsz*i);
      data[i] = vd[i]->vptr;
    }
    H5Tcommit(f, tpath, h5t);
  }
 
  DEBUG("trying to write_uniform_cartesian3Ds(%d, %f, %s, %p, %p, %p, %s, %d, %p, %p, %d)",
        f, t, "Cartesian", orig, space,
        dim, name, h5t, (const void **)data,
        0, H5P_DEFAULT);
  F5Fwrite_uniform_cartesian3Ds(f, t, "Cartesian", orig, space,
                                dim, name, h5t, (const void **)data,
                                0, H5P_DEFAULT);
}
#endif
 
void free_memory(vdata **vd, int n)
{
  int i;
 
  for (i = 0; i < n; i++)
  {
    vinfo *vi = vd[i]->ncref;
 
    free(vd[i]->vptr);
    free(vi->name);
    if (vi->dim) free(vi->dim);
    if (vi->ivar) free(vi->ivar);
  }
  free(vd);
}
 
int main(int argc, char **argv)
{
  int c, nnames[MAXGROUPS], ngroups = 0;
  char **vnames[MAXGROUPS], *vtags[MAXGROUPS];
  char *ofile = "out.f5", *timevar = "time";
  hid_t h5f;
  double t = 0.0, dist;
 
  while ((c = getopt(argc, argv, "v:o:n:t:r:h")) != -1)
  {
    switch (c)
    {
    case 'v':
      if (ngroups >= MAXGROUPS)
        error("exceeded max. number of variable groups");
      vnames[ngroups] = optvars(optarg, nnames+ngroups);
      vtags[ngroups++] = 0;
      break;
    case 'n':
      if (ngroups <= 0)
        error("variable name precedes variable specifier");
      vtags[ngroups-1] = optarg;
      break;
    case 't':
      timevar = optarg;
      break;
    case 'r':
      dist = strtod(optarg, 0);
      break;
    case 'o':
      ofile = optarg;
      break;
    case 'h':
      usage(argv[0], 0);
    default:
      usage(argv[0], 1);
    }
  }
 
  if (optind >= argc) usage(argv[0], 1);
 
  if ((h5f = H5Fcreate(ofile, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)) < 0)
    error("cannot create output F5 file");
 
  while (optind < argc)
  {
    int g, st, ncid;
    vdata **vtab;
    vinfo *uvar, *tinfo;
    char *ifn;
 
    ifn = argv[optind++];
    if ((st = nc_open(ifn, NC_NOWRITE, &ncid) != NC_NOERR))
      error("cannot open NETCDF file");
    tinfo = get_time(ncid, timevar);
 
    for (g = 0; g < ngroups; g++)
    {
      F5_vec3_point_t orig;
      F5_vec3_float_t spc;
 
      if (!(vtab = get_ncvars(ncid, nnames[g], vnames[g])))
        error("failed to get variable info from NETCDF file");
      if (!(uvar = unify(tinfo, vtab, nnames[g])))
        error("failed to unify dimensionalities");
      read_ncvars(ncid, uvar, vtab, nnames[g]);
      dimopt(uvar);
 
      /*
      if (writef5(h5f, t, vtags[g], uvar, vtab, nnames[g]) != F5_SUCCESS)
        error("F5 write failed");
      */
 
      if (!vtags[g]) vtags[g] = "unspecified";
 
      switch (uvar->ndims)
      {
      case 3:
        orig.x = orig.y = orig.z = 0.0;
        spc.x = spc.y = spc.z = dist;
        /*
        if (writef5_cart3d(h5f, t, vtags[g], &orig, &spc,
                           uvar, vtab, nnames[g]) != F5_SUCCESS)
          error("failed to write field on 3D grid");
        */
        switch (nnames[g])
        {
        case 3:
          writef5_uni3d_f3(h5f, t, vtags[g], &orig, &spc, uvar, vtab);
          break;
        case 1:
          writef5_uni3d_scalar(h5f, t, vtags[g], &orig, &spc,
                               uvar, vtab[0]);
          break;
        default:
          error("only scalars or 3D vectors supported");
        }
        break;
      case 1:
        error("1D grids not suppoted yet");
        break;
      default:
        error("unsupported grid rank");
      }
    }
    free_memory(vtab, nnames[g]);
    t += 1;
  }
  F5Xclose(h5f);
 
  return 0;
}