diff -r 000000000000 -r 3c02cce30be8 h1_summary.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/h1_summary.c Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,1525 @@ +#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 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 *mcsec_name = "mcsec"; +static char *history_name = "history"; +static char *nsamples_name = "nsamples"; +static char *hour_name = "hour", *hour_long_name = "hour of day", + *hour_units = "hour"; +static char *cell_method_name = "cell_method"; +/* 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 int mean_hour_dimid; /* special notes */ +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 int stddev_hour_dimid; /* special notes */ +static struct dim *stddev_dim_head = NULL; +static struct var *stddev_var_head = NULL; + +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, int *out_hour_dimidp) +{ + /* + * 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" and "hour" dimensions. + */ + 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; + + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid)); + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY); + ++ndims; + /* Track hour dimid for out files */ + *out_hour_dimidp = dimid; + } + + 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(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid) +{ + /* + * Translate between input dimension IDs and those for the output file. + * For normal time-based variables, add a new dimension for hour of + * day. For metadata variables, do not add new dimensions, even if + * it is time-based. Return (possibly new) number of dimensions. + */ + int i, ndims; + + for (i = ndims = 0; i < in_ndims; i++, ndims++) { + out_dimids[ndims] = in2out[in_dimids[i]]; + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name) + && !metadata) { + ndims++; + out_dimids[ndims] = hour_dimid; + } + } + + 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) +{ + /* + * 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 either to "time: mean" or + * "time: standard_deviation", depending on the output file type. + */ + + int i, natts; + char name[NC_MAX_NAME+1], cmflag = 0; + char *cell_methods[] = { "time: mean", "time: standard_deviation" }; + + 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 + 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, int hour_dimid, struct var **out_var_headp, + 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(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid); + 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); + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + /* + * Just after the "time" variable, add the new "hour" + * variable for hour of day. + */ + if (!strcmp(vnode->name, time_name)) { + ndims = 1; + dimids[0] = hour_dimid; + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid)); + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name)); + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units)); + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0); + ++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 copy_metadata(int in_ncid, struct var *in_var_head, + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head) +{ + int i, j; + float hr[HOURS_PER_DAY]; + struct var *in_vnode, *out_vnode; + void *val; + + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i; + + if (in_var_head) { + in_vnode = in_var_head; + /* + * March through input variables to stuff metadata values into + * the new output files. NOTE: Time-based metadata variables + * should contain only the last (time-slice) value (from all + * input files). + */ + for (j = 0; j == 0 || in_vnode != in_var_head; j++) { + if (in_vnode->metadata) { + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name); + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + /* time is not the first dimension */ +#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); + } + else { + /* time is the first dimension */ +#ifdef DEBUG + printf("Copying last value of \ +time-based metadata variable %s\n", in_vnode->name); +#endif /* DEBUG */ + val = read_timeslice(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, ((input_dim_idx[in_vnode->dimids[0]])->len - 1)); + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0); + } + free(val); + /* + * Just after the "time" variable, write out + * the "hour" variable values. + */ + if (!strcmp(in_vnode->name, time_name)) { + out_vnode = varlist_find_by_name(out_var_head, hour_name); + write_var(out_ncid, out_vnode->ncvarid, + out_vnode->nctype, hr); + } + } + else { +#ifdef DEBUG + printf("Skipping variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + } + in_vnode = in_vnode->next; + } + } + + return; +} + +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples) +{ + char *mean_history_gatt, *stddev_history_gatt, + *mean_prefix="Statistical means from history files:", + *stddev_prefix="Statistical standard deviations from history files:"; + + /* + * 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); + + /* Open input file */ + wrap_nc(nc_open(input_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, &mean_hour_dimid); + stddev_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid, + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid); + /* Define variables and their attributes */ + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head, + MEAN_TYPE); + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head, + STDDEV_TYPE); + /* 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); + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid, + stddev_var_head); + + 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] = 0; + count[0] = 1; /* time */ + start[1] = hr; + count[1] = 1; /* hour */ + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 5) { + start[2] = level; + count[2] = 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, size_t nsamples) +{ + 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[(nifnames-1)], mfname, sfname, flist, nsamples); + + compute_stats(nifnames, ifnames, nsamples); + + /* Close the two output files */ + wrap_nc(nc_close(mean_ncid)); + wrap_nc(nc_close(stddev_ncid)); + + return 0; +}