diff -r d3122367777b -r dd8e6719647b h1_summary_cb.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/h1_summary_cb.c Wed Oct 10 11:59:02 2007 -0400 @@ -0,0 +1,1737 @@ +#include +#include +#include +#include +#include +#include "netcdf.h" + +/* + * Loop through one month of hourly history tapes from the Community Land Model + * and generate monthly statistics (means and standard deviations) of fields + * by hour of day. + */ + +#define MEAN_TYPE 0 +#define STDDEV_TYPE 1 +#define NUM_TYPES 2 + +#define HOURS_PER_DAY 24 +#define MIN_PER_HOUR 60 +#define SEC_PER_MIN 60 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN) + +static char *metadata_vars[] = { + "lon", + "lat", + "lonrof", + "latrof", + "time", + "hour", /* new metadata variable to be added to output files */ + "levsoi", + "levlak", + "edgen", + "edgee", + "edges", + "edgew", + "longxy", + "latixy", + "area", + "areaupsc", + "topo", + "topodnsc", + "landfrac", + "landmask", + "pftmask", + "indxupsc", + "mcdate", + "mcsec", + "mdcur", + "mscur", + "nstep", + "time_bounds", + "date_written", + "time_written" +}; + +struct text_att { + char *name; + char *value; + struct text_att *next; + struct text_att *prev; +}; + +struct dim { + int dimid; + char *name; + size_t len; + struct dim *next; + struct dim *prev; +}; + +struct var { + int ncvarid; + char *name; + nc_type nctype; + int ndims; + int *dimids; + int natts; + char metadata; + char FillFlag; + float FillValue; + struct var *next; + struct var *prev; +}; + +static char *time_name = "time"; +static char *time_bounds_attname = "bounds"; +static char *time_bounds_name = "time_bounds"; +static char *climatology_bounds_attname = "climatology"; +static char *climatology_bounds_name = "climatology_bounds"; +static char *mcdate_name = "mcdate"; +static char *mcsec_name = "mcsec"; +static char *history_name = "history"; +static char *nsamples_name = "nsamples"; +static char *cell_method_name = "cell_method"; +static char *cell_methods_prefix[] = { "time: mean within days time: mean over days", "time: mean within days times: standard_deviation over days" }; +/* input stuff */ +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars); +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid; +static struct text_att *input_gatt_head = NULL; +static struct dim *input_dim_head = NULL, **input_dim_idx; +static struct var *input_var_head = NULL; +/* translation stuff */ +static int *idid2mdid, *idid2sdid; /* dimension IDs */ +/* output stuff */ +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid; +static struct dim *mean_dim_head = NULL; +static struct var *mean_var_head = NULL; +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid; +static struct dim *stddev_dim_head = NULL; +static struct var *stddev_var_head = NULL; +/* output values */ +static double *times; +static double *tbounds; +static int *mcdate; +static int *mcsec; + +char is_metadata(char *name) +{ + int i; + + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++); + + if (i < nmvars) + return 1; + else + return 0; +} + +#if 0 +struct dim *dimlist_find_by_name(struct dim *head, char *name) +{ + int i; + struct dim *p = NULL; + + if (head) { + node = head; + for (i = 0 ; i == 0 || node != head; i++) { + if (!strcmp(node->name, name)) + p = node; + node = node->next; + } + return p; + } + else + return NULL; +} +#endif + +struct var *varlist_find_by_name(struct var *head, char *name) +{ + int i; + struct var *node; + + if (head) { + node = head; + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++) + node = node->next; + if (i && node == head) + return NULL; + else + return node; + } + else + return NULL; +} + +void gattlist_add_node(struct text_att **headp, char *name, char *value) +{ + struct text_att *head, *node; + + head = *headp; + + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) { + perror("gattlist_add_node"); + exit(2); + } + + if (!(node->name = strdup(name))) { + perror("gattlist_add_node"); + exit(2); + } + if (!(node->value = strdup(value))) { + perror("gattlist_add_node"); + exit(2); + } + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len) +{ + struct dim *head, *node; + + head = *headp; + + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) { + perror("dimlist_add_node"); + exit(2); + } + + node->dimid = dimid; + if (!(node->name = strdup(name))) { + perror("dimlist_add_node"); + exit(2); + } + node->len = len; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return node; +} + +void varlist_add_node(struct var **headp, int ncvarid, char *name, + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag, + float FillValue) +{ + int i; + struct var *head, *node; + + head = *headp; + + if (!(node = (struct var *)malloc(sizeof(struct var)))) { + perror("varlist_add_node"); + exit(3); + } + + node->ncvarid = ncvarid; + if (!(node->name = strdup(name))) { + perror("varlist_add_node"); + exit(3); + } + node->nctype = nctype; + node->ndims = ndims; + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) { + perror("varlist_add_node"); + exit(3); + } + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i]; + node->natts = natts; + node->metadata = is_metadata(name); + node->FillFlag = FillFlag; + node->FillValue = FillValue; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +void varlist_print(struct var *headp) +{ + int i, j; + struct var *node; + + printf("Beginning of Variable List\n"); + + if (headp) { + node = headp; + for (j = 0; j == 0 || node != headp; j++) { + printf("Variable %d (ptr=%ld):\n", j, (long)node); + printf("\tncvarid = %d\n", node->ncvarid); + printf("\tname = %s\n", node->name); + printf("\tnctype = %d\n", node->nctype); + printf("\tndims = %d\n", node->ndims); + printf("\tdimids ="); + for (i = 0; i < node->ndims; i++) + printf(" %d", node->dimids[i]); + printf("\n"); + printf("\tnatts = %d\n", node->natts); + printf("\tmetadata = %d\n", (int)node->metadata); + printf("\tnext = %ld\n", (long)node->next); + printf("\tprev = %ld\n", (long)node->prev); + node = node->next; + } + } + else { + printf("\tList undefined\n"); + } + + printf("End of Variable List\n"); + + return; +} + +void wrap_nc(int status) +{ + if (status != NC_NOERR) { + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); + exit(1); + } + + return; +} + +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len; + struct dim **in_dim_idx; + + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) { + perror("*in_dim_idxp"); + exit(3); + } + in_dim_idx = *in_dim_idxp; + + for (i = 0; i < ndims; i++) { + wrap_nc(nc_inq_dim(ncid, i, name, &len)); + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len); + } + + return; +} +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp) +{ + int i; + char name[NC_MAX_NAME+1], *value; + nc_type xtype; + size_t len; + + for (i = 0; i < ngatts; i++) { + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name)); + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len)); + if (xtype != NC_CHAR) { + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name); + exit(2); + } + if (!(value = (char *)malloc((len+1)*sizeof(char)))) { + perror("get_gatts: value"); + exit(3); + } + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value)); + value[(len+1-1)] = '\0'; + gattlist_add_node(gatt_headp, name, value); + free(value); /* because gattlist_add_node() duplicates it */ + } + + return; +} + +void get_vars(int ncid, int nvars, struct var **var_headp) +{ + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts; + size_t len; + float FillValue; + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag; + nc_type xtype, atype; + + for (i = 0; i < nvars; i++) { + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids, + &natts)); + /* Check for _FillValue */ + FillFlag = 0; + FillValue = 0.; + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len) + == NC_NOERR) { + if (atype == NC_FLOAT && len) { + wrap_nc(nc_get_att_float(ncid, i, + fill_att_name, &FillValue)); + FillFlag = 1; + } + } + varlist_add_node(var_headp, i, name, xtype, ndims, dimids, + natts, FillFlag, FillValue); + } + + return; +} + +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid, + size_t nsamples, int **in2outp, int out_ncid, + struct dim **out_dim_headp, int *out_unlimdimidp) +{ + /* + * Define dimensions on new files and build translation tables between + * dimension IDs. + */ + int j, dimid, ndims, *in2out; + size_t len; + struct dim *dnode; + + ndims = 0; + + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) { + perror("put_dimensions:in2outp"); + exit(3); + } + in2out = *in2outp; + + if (in_dim_head) { + dnode = in_dim_head; + for (j = 0; j == 0 || dnode != in_dim_head; j++) { + if (dnode->dimid == in_unlimdimid) + len = NC_UNLIMITED; + else + len = dnode->len; + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid)); + dimlist_add_node(out_dim_headp, dimid, dnode->name, len); + ++ndims; + in2out[dnode->dimid] = dimid; + if (dnode->dimid == in_unlimdimid) + *out_unlimdimidp = dimid; + /* Just after the "time" dimension, add the new + * "nsamples" dimension. */ + if (!strcmp(dnode->name, time_name)) { + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid)); + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples); + ++ndims; + } + dnode = dnode->next; + } + } + else { + fprintf(stderr, "WARNING: No dimensions defined!\n"); + } + + return ndims; +} + +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history) +{ + /* + * Write out global attributes matching those of the input file. + * Change history attribute to the string passed in as an argument. + */ + int j, hflag = 0, ngatts; + char *value; + struct text_att *anode; + + ngatts = 0; + + if (in_gatt_head) { + anode = in_gatt_head; + for (j = 0; j == 0 || anode != in_gatt_head; j++) { + if (!strcmp(anode->name, history_name)) { + value = out_history; + ++hflag; + } + else + value = anode->value; + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value)); + ++ngatts; + anode = anode->next; + } + /* If no history attribute on input, add one to the output */ + if (!hflag) { + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history)); + ++ngatts; + } + } + else { + fprintf(stderr, "WARNING: No global attributes defined!\n"); + } + + return ngatts; +} + +int translate_dimids(int in_ndims, int *in_dimids, int *in2out, int *out_dimids) +{ + /* + * Translate between input dimension IDs and those for the output file. + */ + int i, ndims; + + for (i = ndims = 0; i < in_ndims; i++, ndims++) + out_dimids[ndims] = in2out[in_dimids[i]]; + + return ndims; +} + +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid, + int in_varid, int out_ncid, int out_varid, char **cell_methods) +{ + /* + * Copy the attributes of the input variable to those of the output + * variable. If the variable is not a metadata variable, ensure that + * the cell_method attribute is set correctly, depending on the output + * file type. + */ + + int i, natts; + char name[NC_MAX_NAME+1], cmflag = 0; + + for (i = natts = 0; i < in_natts; i++, natts++) { + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name)); + if (!strcmp(name, cell_method_name) && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + cmflag = 1; + } + else if (!strcmp(name, time_bounds_attname)) + /* Change time_bounds to climatology_bounds */ + wrap_nc(nc_put_att_text(out_ncid, out_varid, climatology_bounds_attname, strlen(climatology_bounds_name), climatology_bounds_name)); + else + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid)); + } + /* + * If no cell_method attribute was specified for a non-metadata + * variable on the input file, add the appropriate cell_method anyway + * on the output file. + */ + if (!cmflag && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + natts++; + } + + return natts; +} + +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head, + int *in2out, int out_ncid, struct var **out_var_headp, + char **cell_methods, int stat_type) +{ + /* + * Define variables on output file and copy attributes from input file. + * Return number of variables defined. + */ + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts; + struct var *vnode; + + nvars = 0; + + if (in_var_head) { + vnode = in_var_head; + /* + * March through input variables creating (mostly) the same + * variables on the output file. + */ + for (j = 0; j == 0 || vnode != in_var_head; j++) { + ndims = translate_dimids(vnode->ndims, vnode->dimids, in2out, dimids); + /* Create all normal data variables, but only the four + * time-based metadata variables: time, + * climatology_bounds, mcdate, and mcsec. */ + if (!strcmp(vnode->name, time_name)) { + /* Magically promote "time" variable to + * double precision on output file */ + wrap_nc(nc_def_var(out_ncid, vnode->name, NC_DOUBLE, ndims, dimids, &varid)); + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods); + varlist_add_node(out_var_headp, varid, vnode->name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + /* Create "climatology_bounds" right after + * time and instead of "time_bounds", with no + * attributes */ + wrap_nc(nc_def_var(out_ncid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, &varid)); + natts = 0; + varlist_add_node(out_var_headp, varid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + } + else if (!vnode->metadata || !strcmp(vnode->name, mcdate_name) || !strcmp(vnode->name, mcsec_name) || (vnode->metadata && strcmp((in_dim_idx[vnode->dimids[0]])->name, time_name))) { + /* If it is not a metadata variable OR it is + * mcdate OR it is mcsec OR (it is a metadata + * variable that does not have time as its + * first dimension), then create it */ + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid)); + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods); + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + } + vnode = vnode->next; + } + } + else { + fprintf(stderr, "WARNING: No variables defined!\n"); + } + + return nvars; +} + +void *alloc_var(nc_type nctype, size_t len) +{ + void *val; + + switch (nctype) { + case NC_FLOAT: + if (!(val = (float *)malloc((len) * sizeof(float)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_INT: + if (!(val = (int *)malloc((len) * sizeof(int)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_DOUBLE: + if (!(val = (double *)malloc((len) * sizeof(double)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_CHAR: + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) { + perror("alloc_var: val"); + exit(3); + } + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx) +{ + int i; + size_t len = (size_t)1; + void *val; + + /* Figure out the total size */ + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len; + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_get_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice) +{ + int i; + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + void *val; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + /* + * Figure out the total size for the slice (start at second dimension) + * and build start/count vectors. + */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + len *= (dim_idx[dimids[i]])->len; + } + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void write_var(int ncid, int varid, nc_type nctype, void *val) +{ + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_put_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice) +{ + int i; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + + /* Build start/count vectors */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + } + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void write_time_var(int ncid, int varid, nc_type nctype, void *time_var) +{ + size_t start[1], count[1]; + + start[0] = 0; + count[0] = 24; + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_vara_float(ncid, varid, start, count, time_var)); + break; + case NC_INT: + wrap_nc(nc_put_vara_int(ncid, varid, start, count, time_var)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, time_var)); + break; + case NC_CHAR: + wrap_nc(nc_put_vara_text(ncid, varid, start, count, time_var)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void write_tbounds(int ncid, int varid, double *tbounds) +{ + size_t start[2], count[2]; + + start[0] = start[1] = 0; + count[0] = 24; count[1] = 2; + + wrap_nc(nc_put_vara_double(ncid, varid, start, count, tbounds)); + + return; +} + +void copy_metadata(int in_ncid, struct var *in_var_head, + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head, + double *times, double *tbounds, int *mcdate, int *mcsec) +{ + int i; + struct var *in_vnode, *out_vnode; + void *val; + + if (!in_var_head) return; + + in_vnode = in_var_head; + /* + * March through input variables to stuff metadata values into + * the new output files. NOTE: The only time-based metadata + * variables, explicitly handled, are: time, climatology_bounds, + * mcdate, and mcsec. These are passed in. + */ + for (i = 0; i == 0 || in_vnode != in_var_head; i++) { + if (in_vnode->metadata) { + /* Find the corresponding output variable */ + if (!strcmp(in_vnode->name, time_bounds_name)) + out_vnode = varlist_find_by_name(out_var_head, climatology_bounds_name); + else + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name); + if (out_vnode) { + /* Read values from input and write to output */ + if (strcmp((in_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + /* "time" is not the first dimension, + * so just copy the variable values */ +#ifdef DEBUG + printf("Copying metadata variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx); + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val); + free(val); + } + else { + /* "time" is the first dimension, so + * explicitly handle the "time", + * "climatology_bounds", "mcdate", and + * "mcsec" while dropping all others */ + if (!strcmp(out_vnode->name, time_name)) { + /* Force time stamps to be + * those computed previously, + * and remember these are now + * double precision. */ +#ifdef DEBUG + printf("Writing metadata variable %s\n", out_vnode->name); +#endif /* DEBUG */ + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, times); + } + if (!strcmp(out_vnode->name, climatology_bounds_name)) { + /* time_bounds are now + * climatology_bounds and have + * pre-computed values */ +#ifdef DEBUG + printf("Writing metadata variable %s\n", out_vnode->name); +#endif /* DEBUG */ + write_tbounds(out_ncid, out_vnode->ncvarid, tbounds); + } + else if (!strcmp(out_vnode->name, mcdate_name)) { + /* mcdate is pre-computed */ +#ifdef DEBUG + printf("Writing metadata variable %s\n", out_vnode->name); +#endif /* DEBUG */ + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcdate); + } + else if (!strcmp(out_vnode->name, mcsec_name)) { + /* mcsec is pre-computed */ +#ifdef DEBUG + printf("Writing metadata variable %s\n", out_vnode->name); +#endif /* DEBUG */ + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcsec); + } + } + } + else { +#ifdef DEBUG + printf("Dropping metadata variable %s\n", in_vnode->name); +#endif /* DEBUG */ + } + } + else { +#ifdef DEBUG + printf("Skipping non-metadata variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + } +#ifdef DEBUG + printf ("Done\n"); +#endif /* DEBUG */ + in_vnode = in_vnode->next; + } + + return; +} + +void get_time_bounds(char *first_fname, char *last_fname, double *times, double *tbounds, int *mcdate, int *mcsec) +{ + int i, time_dimid, time_bounds_varid, mcdate_varid, mcsec_varid; + size_t time_len, time_bounds_start[2], time_bounds_count[2], + mc_start[1], mc_count[1]; + double bnd1[24], bnd2[24], day1, day2; + + /* Open first input file */ + wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid)); + /* Get dimension ID of the time dimension */ + wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid)); + /* Get length of time dimension */ + wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len)); + /* Get variable ID of time_bounds variable */ + wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid)); + /* Set start/count for reading the left values of time_bounds from the + * first 24 timeslices of the first file */ + time_bounds_start[0] = time_bounds_start[1] = 0; + time_bounds_count[0] = 24; time_bounds_count[1] = 1; + /* Read the values */ + wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid, + time_bounds_start, time_bounds_count, bnd1)); + /* If the first and last file are not the same, close the first one and + * open the second one */ + if (strcmp(first_fname, last_fname)) { + /* Close the first input file */ + wrap_nc(nc_close(input_ncid)); + /* Open the last input file */ + wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid)); + /* Get dimension ID of the time dimension */ + wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid)); + /* Get length of time dimension */ + wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len)); + /* Get variable ID of time_bounds variable */ + wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, + &time_bounds_varid)); + } + /* Set start/count for reading the right values of time_bounds from the + * last 24 timeslices of the last file */ + time_bounds_start[0] = time_len - 24; time_bounds_start[1] = 1; + time_bounds_count[0] = 24; time_bounds_count[1] = 1; + /* Read the value */ + wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid, + time_bounds_start, time_bounds_count, bnd2)); + /* Read the last 24 mcdate and mcsec values */ + wrap_nc(nc_inq_varid(input_ncid, mcdate_name, &mcdate_varid)); + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &mcsec_varid)); + mc_start[0] = time_len - 24; + mc_count[0] = 24; + wrap_nc(nc_get_vara_int(input_ncid, mcdate_varid, mc_start, mc_count, + mcdate)); + wrap_nc(nc_get_vara_int(input_ncid, mcsec_varid, mc_start, mc_count, + mcsec)); + + /* Close the last input file */ + wrap_nc(nc_close(input_ncid)); + + /* Store time bounds and compute new time stamps as midpoints */ + for (i = 0; i < 24; i++) { + tbounds[(i*2)] = bnd1[i]; + tbounds[(i*2+1)] = bnd2[i]; + /* Pick one (integer) day near the midpoint of the bounds and + * add in the fraction from the right bounds, so for January + * the first hour-of-day value would be 15.0 (plus some year + * offset) */ + day1 = floor(bnd1[i]); + day2 = floor(bnd2[i]); + times[i] = floor((day1 + day2) / 2.0) + (bnd2[i]-day2); + } + + return; +} + +char **make_cell_methods(size_t nsamples) +{ + int i; + size_t len; + char **cell_methods; + + if (!(cell_methods = (char **)malloc((NUM_TYPES * sizeof(char *))))) { + perror("make_cell_methods:cell_methods"); + exit(2); + } + for (i = 0; i < NUM_TYPES; i++) { + len = strlen(cell_methods_prefix[i]) + 48; + if (!(cell_methods[i] = (char *)malloc((len * sizeof(char))))) { + perror("make_cell_methods:cell_methods[i]"); + exit(2); + } + sprintf(cell_methods[i], "%s (comment: %ld samples)", + cell_methods_prefix[i], (long)nsamples); + } + return cell_methods; +} + +void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples) +{ + int i; + char *mean_history_gatt, *stddev_history_gatt, + *mean_prefix="Statistical means from history files:", + *stddev_prefix="Statistical standard deviations from history files:", + **cell_methods; + + /* + * Construct strings for history global attributes for the two output + * files. + */ + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:mean_history_gatt"); + exit(2); + } + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:stddev_history_gatt"); + exit(2); + } + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist); + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist); + + /* Allocate space for time values for each hour of the day */ + if (!(times = (double *)malloc((24 * sizeof(double))))) { + perror("open_inout:times"); + exit(2); + } + /* Allocate space for two time bounds values for each hour of the day */ + if (!(tbounds = (double *)malloc((24 * 2 * sizeof(double))))) { + perror("open_inout:tbounds"); + exit(2); + } + /* Allocate space for current date and seconds for each hour of day */ + if (!(mcdate = (int *)malloc((24 * sizeof(int))))) { + perror("open_inout:mcdate"); + exit(2); + } + if (!(mcsec = (int *)malloc((24 * sizeof(int))))) { + perror("open_inout:mcsec"); + exit(2); + } + + + /* + * Explicitly handle "time", "time_bounds" -> "climatology_bounds", + * mcdate, and mcsec metadata variables by constructing their values + * before doing anything else. + */ + get_time_bounds(first_fname, last_fname, times, tbounds, mcdate, mcsec); +#ifdef DEBUG + for (i = 0; i < 24; i++) + fprintf(stderr, "Hour of day=%d, time=%lf, tbounds=%lf,%lf mcdate=%d mcsec=%d\n", i, times[i], tbounds[(i*2)], tbounds[(i*2+1)], mcdate[i], mcsec[i]); +#endif /* DEBUG */ + + /* Open last input file */ + wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid)); + /* Inquire about number of dimensions, variables, global attributes + * and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts, + &input_unlimdimid)); + /* Grab dimension IDs and lengths from input file */ + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx); + /* Grab desired global attributes from input file */ + get_gatts(input_ncid, input_ngatts, &input_gatt_head); + /* Grab variable IDs and attributes from input file */ + get_vars(input_ncid, input_nvars, &input_var_head); +#ifdef DEBUG + varlist_print(input_var_head); +#endif /* DEBUG */ + /* Create netCDF files for monthly means and standard deviations */ + /* Create new netCDF files */ + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid)); + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid)); + /* Define dimensions */ + mean_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2mdid, mean_ncid, + &mean_dim_head, &mean_unlimdimid); + stddev_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid, + &stddev_dim_head, &stddev_unlimdimid); + /* Define cell_methods[], including number of samples */ + cell_methods = make_cell_methods(nsamples); + /* Define variables and their attributes */ + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2mdid, mean_ncid, &mean_var_head, cell_methods, + MEAN_TYPE); + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2sdid, stddev_ncid, &stddev_var_head, cell_methods, + STDDEV_TYPE); +#ifdef DEBUG + printf("Mean variables:\n"); + varlist_print(mean_var_head); + printf("Stddev variables:\n"); + varlist_print(stddev_var_head); +#endif /* DEBUG */ + /* Store global attributes */ + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt); + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid, + stddev_history_gatt); + /* End define mode */ + wrap_nc(nc_enddef(mean_ncid)); + wrap_nc(nc_enddef(stddev_ncid)); + /* Write out metdata variables that do not depend on "time" */ + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid, + mean_var_head, times, tbounds, mcdate, mcsec); + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid, + stddev_var_head, times, tbounds, mcdate, mcsec); + + wrap_nc(nc_close(input_ncid)); + + return; +} + +size_t count_samples(int nifnames, char **input_fnames) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len, cnt = 0; + + for (i = 0; i < nifnames; i++) { + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } +#ifdef DEBUG + printf("%ld samples in %s\n", (long)len, input_fnames[i]); +#endif /* DEBUG */ + wrap_nc(nc_close(input_ncid)); + cnt += len; + } + return cnt; +} + +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp) +{ + /* + * Allocate space for arrays of means and standard deviations by + * hour of day. + */ + int i; + size_t **cell_samples; + double **means, **stddevs; + + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*meansp"); + exit(2); + } + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*stddevsp"); + exit(2); + } + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) { + perror("alloc_means_stddevs:*cell_samplesp"); + exit(2); + } + + means = *meansp; + stddevs = *stddevsp; + cell_samples = *cell_samplesp; + + for (i = 0; i < HOURS_PER_DAY; i++) { + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:means[i]"); + exit(2); + } + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:stddevs[i]"); + exit(2); + } + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) { + perror("alloc_means_stddevs:cell_samples[i]"); + exit(2); + } + } + + return; +} + +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue) +{ + int hr, i, j, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + means[hr][pos] = FillValue; + stddevs[hr][pos] = FillValue; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples) +{ + /* + * Free space from arrays of means and standard deviations by + * hour of day. + */ + int i; + + for (i = 0; i < HOURS_PER_DAY; i++) { + free(means[i]); + free(stddevs[i]); + free(cell_samples[i]); + } + + free(means); + free(stddevs); + free(cell_samples); + + return; +} + +void add_to_means_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = (double)val[pos]; + else + means[hr][pos] += (double)val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_means_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = val[pos]; + else + means[hr][pos] += val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + means[hr][pos] /= (double)cell_samples[hr][pos]; + } + } + } + } + + return; +} + +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - (double)val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + stddevs[hr][pos] /= (double)cell_samples[hr][pos]; + stddevs[hr][pos] = sqrt(stddevs[hr][pos]); + } + } + } + } + + return; +} + +void read_and_compute(int nifnames, char **input_fnames, size_t d1, size_t d2, char *var_name, nc_type nctype, int level, int ndims, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples) +{ + int i, ts; + int varid, *mcsec; + size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + char name[NC_MAX_NAME+1]; + void *val; + + /* Allocate space for timeslice */ + val = alloc_var(nctype, (d1 * d2)); + + for (i = 0; i < nifnames; i++) { +#ifdef DEBUG + printf("\tOpening %s", input_fnames[i]); + if (ndims > 3) printf(" to retrieve level %d", level); + if (stat_type == MEAN_TYPE) + printf(" for computing means\n"); + else + printf(" for computing stddevs\n"); +#endif /* DEBUG */ + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } + /* Get seconds of day variable */ + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid)); + if (!(mcsec = (int *)malloc(len * sizeof(int)))) { + perror("read_and_compute:mcsec"); + exit(2); + } + wrap_nc(nc_get_var_int(input_ncid, varid, mcsec)); + /* Get variable ID for requested variable */ + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid)); + /* Retrieve time slice of desired variable */ + for (ts = 0; ts < len; ts++) { + start[0] = ts; + count[0] = 1; + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 4) { + start[1] = level; + count[1] = 1; + } + switch(nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val)); + if (stat_type == MEAN_TYPE) + add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val)); + if (stat_type == MEAN_TYPE) + add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + /* Free mcsec vector */ + free(mcsec); + /* Close input file */ + wrap_nc(nc_close(input_ncid)); + } + /* Divide sums by number of samples */ + if (stat_type == MEAN_TYPE) + divide_means(d1, d2, means, cell_samples); + else + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples); + + /* Free working variable array */ + free(val); + + return; +} + +float *double_to_float(size_t len, double *dvar) +{ + int i; + float *fvar; + + if (!(fvar = (float *)malloc(len * sizeof(float)))) { + perror("double_to_float:fvar"); + exit(2); + } + + for (i = 0; i < len; i++) + fvar[i] = (float)dvar[i]; + + return fvar; +} + +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2, + int level, int ndims, double **var_hours) +{ + /* Output dimensions should be one larger than input dimensions */ + int i, hr; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + float *fvar = NULL; + + if (nctype == NC_FLOAT) { + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) { + perror("write_var_hours:fvar"); + exit(2); + } + } + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + start[0] = hr; + count[0] = 1; /* time */ + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 4) { + start[1] = level; + count[1] = 1; + } + switch (nctype) { + case NC_FLOAT: + for (i = 0; i < (d1 * d2); i++) + fvar[i] = (float)var_hours[hr][i]; + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr])); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + if (nctype == NC_FLOAT) + free(fvar); + + return; +} + +void compute_stats(int nifnames, char **input_fnames) +{ + int j, k, nlevels; + size_t d1, d2, **cell_samples; + double **means, **stddevs; + struct var *in_vnode, *out_vnode; + + if (input_var_head) { + in_vnode = input_var_head; + /* March through non-metadata variables to compute stats */ + for (j = 0; j == 0 || in_vnode != input_var_head; j++) { + if (!in_vnode->metadata) { + /* Make sure time is really the first dimension */ + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name); + exit(3); + } + /* Figure out input dimensions */ + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) { + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims); + exit(3); + } + /* Assume 2D output; find dimensions */ + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len; + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len; + if (in_vnode->ndims == 3) + nlevels = 1; + else + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len; + /* Allocate working space for means and stddevs */ + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples); + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue); + printf("Computing statistics for %s\n", + in_vnode->name); + /* Compute means and stddevs, then write them */ + for (k = 0; k < nlevels; k++) { + /* Read and compute means */ + read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples); + /* Find corresponding output variable on the mean output file */ + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name); + /* Write out the means for this variable */ + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means); + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */ + reset_cell_samples(d1, d2, cell_samples); + /* Read and compute stddevs using means */ + read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples); + /* Find corresponding output variable on the stddev output file */ + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name); + /* Write out stddevs for this variable */ + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs); + } + + /* Free working space for means and stddevs */ + free_means_stddevs(means, stddevs, cell_samples); + } + in_vnode = in_vnode->next; + } + } + return; +} + +void usage(char *program) +{ + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program); + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n"); + exit(4); +} + +int main(int argc, char **argv) +{ + int i, ic, nifnames; + size_t len, pos, nsamples; + char *mfname, *sfname, **ifnames, *flist; + + ifnames = NULL; + mfname = sfname = flist = NULL; + nifnames = 0; + +#ifdef DEBUG + printf("Number of metadata variables (nmvars) = %d\n", nmvars); + fflush(stdout); +#endif /* DEBUG */ + + /* check command-line flags and switches */ + while ((ic = getopt(argc, argv, "m:s:")) != -1) { + switch(ic) { + case 'm': + if (!(mfname = strdup(optarg))) { + perror("mfname"); + exit(2); + } + break; + case 's': + if (!(sfname = strdup(optarg))) { + perror("sfname"); + exit(2); + } + break; + } + } + + if (!mfname) { + fprintf(stderr, "Output file name for writing means required.\n"); + usage(argv[0]); + } + if (!sfname) { + fprintf(stderr, "Output file name for writing standard deviations required.\n"); + usage(argv[0]); + } + + if (optind < argc) { + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) { + perror("ifnames"); + exit(2); + } + for (i = optind; i < argc; i++) { + if (!(ifnames[nifnames++] = strdup(argv[i]))) { + perror("ifnames[nifnames++]"); + exit(2); + } + } + } + else { + fprintf(stderr, "At least one input file name is required.\n"); + usage(argv[0]); + } + + + /* + * Build list of input files to be included in the global history + * attribute on the two outputs files. + */ + for (i = len = 0; i < nifnames; i++) + len += strlen(ifnames[i]); + len += nifnames + 1; /* space in front of every name + null terminator */ + if (!(flist = (char *)malloc(len * sizeof(char)))) { + perror("flist"); + exit(2); + } + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++) + sprintf(&flist[pos], " %s", ifnames[i]); +#ifdef DEBUG + printf("flist=%s\n", flist); + fflush(stdout); +#endif /* DEBUG */ + + /* Open every file to count up the number of time samples. */ + nsamples = count_samples(nifnames, ifnames); +#ifdef DEBUG + printf("Number of samples across all files = %ld\n", (long)nsamples); + fflush(stdout); +#endif /* DEBUG */ + + /* + * Open the *last* input file and the two output files (for mean and + * standard deviation). Define dimensions, variables, and attributes + * for the two output files based on the *last* input files. The + * last file is used because some metadata variables will contain + * only values from the last time slice from the period over which + * the statistics are computed. + */ + open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist, + nsamples); + + compute_stats(nifnames, ifnames); + + /* Close the two output files */ + wrap_nc(nc_close(mean_ncid)); + wrap_nc(nc_close(stddev_ncid)); + + return 0; +}