Fixed h1_summary and h1_summary2 to correctly construct time_bounds values.
h1_summary and h1_summary2 previously used the time_bounds from the last time
value from the last input file when summarizing, suggesting that the field
values were appropriate only over that short time range instead of the complete
time period over which the statistics were calculated. C-LAMP Experiment 1
runs used the previous code, so the time_bounds were incorrect in the
statistical summaries produced. C-LAMP Experiment 2 runs will use this new
code for production of statistical summaries.
9 * Loop through one month of hourly history tapes from the Community Land Model
10 * and generate monthly statistics (means and standard deviations) of fields
17 #define HOURS_PER_DAY 24
18 #define MIN_PER_HOUR 60
19 #define SEC_PER_MIN 60
20 #define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
22 static char *metadata_vars[] = {
28 "hour", /* new metadata variable to be added to output files */
58 struct text_att *next;
59 struct text_att *prev;
84 static char *time_name = "time";
85 static char *time_bounds_name = "time_bounds";
86 static char *mcsec_name = "mcsec";
87 static char *history_name = "history";
88 static char *nsamples_name = "nsamples";
89 static char *hour_name = "hour", *hour_long_name = "hour of day",
91 static char *cell_method_name = "cell_method";
93 static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
94 static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
95 static struct text_att *input_gatt_head = NULL;
96 static struct dim *input_dim_head = NULL, **input_dim_idx;
97 static struct var *input_var_head = NULL;
98 /* translation stuff */
99 static int *idid2mdid, *idid2sdid; /* dimension IDs */
101 static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
102 static int mean_hour_dimid; /* special notes */
103 static struct dim *mean_dim_head = NULL;
104 static struct var *mean_var_head = NULL;
105 static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
106 static int stddev_hour_dimid; /* special notes */
107 static struct dim *stddev_dim_head = NULL;
108 static struct var *stddev_var_head = NULL;
110 char is_metadata(char *name)
114 for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
123 struct dim *dimlist_find_by_name(struct dim *head, char *name)
126 struct dim *p = NULL;
130 for (i = 0 ; i == 0 || node != head; i++) {
131 if (!strcmp(node->name, name))
142 struct var *varlist_find_by_name(struct var *head, char *name)
149 for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
151 if (i && node == head)
160 void gattlist_add_node(struct text_att **headp, char *name, char *value)
162 struct text_att *head, *node;
166 if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
167 perror("gattlist_add_node");
171 if (!(node->name = strdup(name))) {
172 perror("gattlist_add_node");
175 if (!(node->value = strdup(value))) {
176 perror("gattlist_add_node");
188 node->prev = head->prev;
189 /* set this after setting node->prev from here */
191 /* set this after having set node->prev */
192 node->prev->next = node;
198 struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
200 struct dim *head, *node;
204 if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
205 perror("dimlist_add_node");
210 if (!(node->name = strdup(name))) {
211 perror("dimlist_add_node");
224 node->prev = head->prev;
225 /* set this after setting node->prev from here */
227 /* set this after having set node->prev */
228 node->prev->next = node;
234 void varlist_add_node(struct var **headp, int ncvarid, char *name,
235 nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
239 struct var *head, *node;
243 if (!(node = (struct var *)malloc(sizeof(struct var)))) {
244 perror("varlist_add_node");
248 node->ncvarid = ncvarid;
249 if (!(node->name = strdup(name))) {
250 perror("varlist_add_node");
253 node->nctype = nctype;
255 if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
256 perror("varlist_add_node");
259 for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
261 node->metadata = is_metadata(name);
262 node->FillFlag = FillFlag;
263 node->FillValue = FillValue;
273 node->prev = head->prev;
274 /* set this after setting node->prev from here */
276 /* set this after having set node->prev */
277 node->prev->next = node;
283 void varlist_print(struct var *headp)
288 printf("Beginning of Variable List\n");
292 for (j = 0; j == 0 || node != headp; j++) {
293 printf("Variable %d (ptr=%ld):\n", j, (long)node);
294 printf("\tncvarid = %d\n", node->ncvarid);
295 printf("\tname = %s\n", node->name);
296 printf("\tnctype = %d\n", node->nctype);
297 printf("\tndims = %d\n", node->ndims);
298 printf("\tdimids =");
299 for (i = 0; i < node->ndims; i++)
300 printf(" %d", node->dimids[i]);
302 printf("\tnatts = %d\n", node->natts);
303 printf("\tmetadata = %d\n", (int)node->metadata);
304 printf("\tnext = %ld\n", (long)node->next);
305 printf("\tprev = %ld\n", (long)node->prev);
310 printf("\tList undefined\n");
313 printf("End of Variable List\n");
318 void wrap_nc(int status)
320 if (status != NC_NOERR) {
321 fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
328 void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
331 char name[NC_MAX_NAME+1];
333 struct dim **in_dim_idx;
335 if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
336 perror("*in_dim_idxp");
339 in_dim_idx = *in_dim_idxp;
341 for (i = 0; i < ndims; i++) {
342 wrap_nc(nc_inq_dim(ncid, i, name, &len));
343 in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
348 void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
351 char name[NC_MAX_NAME+1], *value;
355 for (i = 0; i < ngatts; i++) {
356 wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
357 wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
358 if (xtype != NC_CHAR) {
359 fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
362 if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
363 perror("get_gatts: value");
366 wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
367 value[(len+1-1)] = '\0';
368 gattlist_add_node(gatt_headp, name, value);
369 free(value); /* because gattlist_add_node() duplicates it */
375 void get_vars(int ncid, int nvars, struct var **var_headp)
377 int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
380 char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
381 nc_type xtype, atype;
383 for (i = 0; i < nvars; i++) {
384 wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
386 /* Check for _FillValue */
389 if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
391 if (atype == NC_FLOAT && len) {
392 wrap_nc(nc_get_att_float(ncid, i,
393 fill_att_name, &FillValue));
397 varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
398 natts, FillFlag, FillValue);
404 int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
405 size_t nsamples, int **in2outp, int out_ncid,
406 struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)
409 * Define dimensions on new files and build translation tables between
412 int j, dimid, ndims, *in2out;
418 if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
419 perror("put_dimensions: in2outp");
426 for (j = 0; j == 0 || dnode != in_dim_head; j++) {
427 if (dnode->dimid == in_unlimdimid)
431 wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
432 dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
434 in2out[dnode->dimid] = dimid;
435 if (dnode->dimid == in_unlimdimid)
436 *out_unlimdimidp = dimid;
438 * Just after the "time" dimension, add the new
439 * "nsamples" and "hour" dimensions.
441 if (!strcmp(dnode->name, time_name)) {
442 wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
443 dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
446 wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
447 dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
449 /* Track hour dimid for out files */
450 *out_hour_dimidp = dimid;
457 fprintf(stderr, "WARNING: No dimensions defined!\n");
463 int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
466 * Write out global attributes matching those of the input file.
467 * Change history attribute to the string passed in as an argument.
469 int j, hflag = 0, ngatts;
471 struct text_att *anode;
476 anode = in_gatt_head;
477 for (j = 0; j == 0 || anode != in_gatt_head; j++) {
478 if (!strcmp(anode->name, history_name)) {
483 value = anode->value;
484 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
488 /* If no history attribute on input, add one to the output */
490 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
495 fprintf(stderr, "WARNING: No global attributes defined!\n");
501 int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)
504 * Translate between input dimension IDs and those for the output file.
505 * For normal time-based variables, add a new dimension for hour of
506 * day. For metadata variables, do not add new dimensions, even if
507 * it is time-based. Return (possibly new) number of dimensions.
511 for (i = ndims = 0; i < in_ndims; i++, ndims++) {
512 out_dimids[ndims] = in2out[in_dimids[i]];
513 if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
516 out_dimids[ndims] = hour_dimid;
523 int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
524 int in_varid, int out_ncid, int out_varid)
527 * Copy the attributes of the input variable to those of the output
528 * variable. If the variable is not a metadata variable, ensure that
529 * the cell_method attribute is set either to "time: mean" or
530 * "time: standard_deviation", depending on the output file type.
534 char name[NC_MAX_NAME+1], cmflag = 0;
535 char *cell_methods[] = { "time: mean", "time: standard_deviation" };
537 for (i = natts = 0; i < in_natts; i++, natts++) {
538 wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
539 if (!strcmp(name, cell_method_name) && !metadata) {
540 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
544 wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
547 * If no cell_method attribute was specified for a non-metadata
548 * variable on the input file, add the appropriate cell_method anyway
549 * on the output file.
551 if (!cmflag && !metadata) {
552 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
559 int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
560 int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp,
564 * Define variables on output file and copy attributes from input file.
565 * Return number of variables defined.
567 int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
575 * March through input variables creating (mostly) the same
576 * variables on the output file.
578 for (j = 0; j == 0 || vnode != in_var_head; j++) {
579 ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
580 wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
581 natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
582 varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
585 * Just after the "time" variable, add the new "hour"
586 * variable for hour of day.
588 if (!strcmp(vnode->name, time_name)) {
590 dimids[0] = hour_dimid;
591 wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
592 wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
593 wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
594 varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
602 fprintf(stderr, "WARNING: No variables defined!\n");
608 void *alloc_var(nc_type nctype, size_t len)
614 if (!(val = (float *)malloc((len) * sizeof(float)))) {
615 perror("alloc_var: val");
620 if (!(val = (int *)malloc((len) * sizeof(int)))) {
621 perror("alloc_var: val");
626 if (!(val = (double *)malloc((len) * sizeof(double)))) {
627 perror("alloc_var: val");
632 if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
633 perror("alloc_var: val");
638 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
645 void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
648 size_t len = (size_t)1;
651 /* Figure out the total size */
652 for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
654 val = alloc_var(nctype, len);
658 wrap_nc(nc_get_var_float(ncid, varid, val));
661 wrap_nc(nc_get_var_int(ncid, varid, val));
664 wrap_nc(nc_get_var_double(ncid, varid, val));
667 wrap_nc(nc_get_var_text(ncid, varid, val));
670 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
677 void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
680 size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
683 /* Make sure time is really the first dimension */
684 if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
685 fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
689 * Figure out the total size for the slice (start at second dimension)
690 * and build start/count vectors.
694 for (i = 1; i < ndims; i++) {
696 count[i] = (dim_idx[dimids[i]])->len;
697 len *= (dim_idx[dimids[i]])->len;
700 val = alloc_var(nctype, len);
704 wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
707 wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
710 wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
713 wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
716 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
723 void write_var(int ncid, int varid, nc_type nctype, void *val)
727 wrap_nc(nc_put_var_float(ncid, varid, val));
730 wrap_nc(nc_put_var_int(ncid, varid, val));
733 wrap_nc(nc_put_var_double(ncid, varid, val));
736 wrap_nc(nc_put_var_text(ncid, varid, val));
739 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
746 void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
749 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
751 /* Make sure time is really the first dimension */
752 if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
753 fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
757 /* Build start/count vectors */
760 for (i = 1; i < ndims; i++) {
762 count[i] = (dim_idx[dimids[i]])->len;
767 wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
770 wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
773 wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
776 wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
779 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
786 void copy_metadata(int in_ncid, struct var *in_var_head,
787 struct dim **in_dim_idx, int out_ncid, struct var *out_var_head,
791 float hr[HOURS_PER_DAY];
792 struct var *in_vnode, *out_vnode;
795 for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
798 in_vnode = in_var_head;
800 * March through input variables to stuff metadata values into
801 * the new output files. NOTE: Time-based metadata variables
802 * should contain only the last (time-slice) value (from all
805 for (j = 0; j == 0 || in_vnode != in_var_head; j++) {
806 if (in_vnode->metadata) {
807 out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
808 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
809 /* time is not the first dimension */
811 printf("Copying metadata variable %s\n",
814 val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
815 write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
818 /* time is the first dimension */
820 printf("Copying last value of \
821 time-based metadata variable %s\n", in_vnode->name);
823 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));
824 if (!strcmp(in_vnode->name, time_bounds_name))
825 write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, tbounds, 0);
827 write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
831 * Just after the "time" variable, write out
832 * the "hour" variable values.
834 if (!strcmp(in_vnode->name, time_name)) {
835 out_vnode = varlist_find_by_name(out_var_head, hour_name);
836 write_var(out_ncid, out_vnode->ncvarid,
837 out_vnode->nctype, hr);
842 printf("Skipping variable %s\n",
846 in_vnode = in_vnode->next;
853 void get_time_bounds(char *first_fname, char *last_fname, double *tbounds)
855 int time_dimid, time_bounds_varid;
856 size_t time_len, time_bounds_index[2];
859 /* Open first input file */
860 wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
861 /* Get dimension ID of the time dimension */
862 wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
863 /* Get length of time dimension */
864 wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
865 /* Get variable ID of time_bounds variable */
866 wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
867 /* Set index for reading the first value of time_bounds from the
868 * first timeslice of the first file */
869 time_bounds_index[0] = time_bounds_index[1] = 0;
871 wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
872 time_bounds_index, &bnd1));
873 /* If the first and last file are not the same, close the first one and
874 * open the second one */
875 if (strcmp(first_fname, last_fname)) {
876 /* Close the first input file */
877 wrap_nc(nc_close(input_ncid));
878 /* Open the last input file */
879 wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
880 /* Get dimension ID of the time dimension */
881 wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
882 /* Get length of time dimension */
883 wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
884 /* Get variable ID of time_bounds variable */
885 wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
886 &time_bounds_varid));
888 /* Set index for reading the second value of time_bounds from the
889 * last timeslice of the last file */
890 time_bounds_index[0] = time_len - 1; time_bounds_index[1] = 1;
892 wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
893 time_bounds_index, &bnd2));
895 /* Close the last input file */
896 wrap_nc(nc_close(input_ncid));
904 void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
906 char *mean_history_gatt, *stddev_history_gatt,
907 *mean_prefix="Statistical means from history files:",
908 *stddev_prefix="Statistical standard deviations from history files:";
912 * Construct strings for history global attributes for the two output
915 if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
916 perror("open_inout:mean_history_gatt");
919 if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
920 perror("open_inout:stddev_history_gatt");
923 sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
924 sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
926 /* The two time_bounds values must be handled explicitly by obtaining
927 * the first one from the first file and the last one from the last
928 * file. These are then passed to copy_metadata() below. */
929 get_time_bounds(first_fname, last_fname, tbounds);
931 fprintf(stderr, "Got back tbounds=%lf,%lf\n", tbounds[0], tbounds[1]);
934 /* Open last input file */
935 wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
936 /* Inquire about number of dimensions, variables, global attributes
937 * and the ID of the unlimited dimension
939 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
941 /* Grab dimension IDs and lengths from input file */
942 get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
943 /* Grab desired global attributes from input file */
944 get_gatts(input_ncid, input_ngatts, &input_gatt_head);
945 /* Grab variable IDs and attributes from input file */
946 get_vars(input_ncid, input_nvars, &input_var_head);
948 varlist_print(input_var_head);
950 /* Create netCDF files for monthly means and standard deviations */
951 /* Create new netCDF files */
952 wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
953 wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
954 /* Define dimensions */
955 mean_ndims = put_dimensions(input_dim_head, input_ndims,
956 input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
957 &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
958 stddev_ndims = put_dimensions(input_dim_head, input_ndims,
959 input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
960 &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
961 /* Define variables and their attributes */
962 mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
963 idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
965 stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
966 idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
968 /* Store global attributes */
969 mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
970 stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
971 stddev_history_gatt);
972 /* End define mode */
973 wrap_nc(nc_enddef(mean_ncid));
974 wrap_nc(nc_enddef(stddev_ncid));
975 /* Write out metdata variables that do not depend on "time" */
976 copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
977 mean_var_head, tbounds);
978 copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
979 stddev_var_head, tbounds);
981 wrap_nc(nc_close(input_ncid));
986 size_t count_samples(int nifnames, char **input_fnames)
989 char name[NC_MAX_NAME+1];
992 for (i = 0; i < nifnames; i++) {
993 /* Open input file */
994 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
996 * Inquire about number of dimensions, variables, global
997 * attributes and the ID of the unlimited dimension
999 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1000 &input_ngatts, &input_unlimdimid));
1001 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1002 if (strcmp(name, time_name)) {
1003 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1007 printf("%ld samples in %s\n", (long)len, input_fnames[i]);
1009 wrap_nc(nc_close(input_ncid));
1015 void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
1018 * Allocate space for arrays of means and standard deviations by
1022 size_t **cell_samples;
1023 double **means, **stddevs;
1025 if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1026 perror("alloc_means_stddevs:*meansp");
1029 if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1030 perror("alloc_means_stddevs:*stddevsp");
1033 if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
1034 perror("alloc_means_stddevs:*cell_samplesp");
1039 stddevs = *stddevsp;
1040 cell_samples = *cell_samplesp;
1042 for (i = 0; i < HOURS_PER_DAY; i++) {
1043 if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1044 perror("alloc_means_stddevs:means[i]");
1047 if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1048 perror("alloc_means_stddevs:stddevs[i]");
1051 if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
1052 perror("alloc_means_stddevs:cell_samples[i]");
1060 void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
1064 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1065 for (i = 0; i < d1; i++) {
1066 #pragma _CRI concurrent
1067 for (j = 0; j < d2; j++) {
1069 means[hr][pos] = FillValue;
1070 stddevs[hr][pos] = FillValue;
1071 cell_samples[hr][pos] = 0;
1079 void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
1083 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1084 for (i = 0; i < d1; i++) {
1085 #pragma _CRI concurrent
1086 for (j = 0; j < d2; j++) {
1088 cell_samples[hr][pos] = 0;
1096 void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
1099 * Free space from arrays of means and standard deviations by
1104 for (i = 0; i < HOURS_PER_DAY; i++) {
1107 free(cell_samples[i]);
1117 void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
1118 char FillFlag, float FillValue, double **means, size_t **cell_samples)
1122 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1124 for (i = 0; i < d1; i++) {
1125 #pragma _CRI concurrent
1126 for (j = 0; j < d2; j++) {
1128 if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1129 if (cell_samples[hr][pos] == 0)
1130 means[hr][pos] = (double)val[pos];
1132 means[hr][pos] += (double)val[pos];
1133 ++cell_samples[hr][pos];
1141 void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
1142 char FillFlag, float FillValue, double **means, size_t **cell_samples)
1146 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1148 for (i = 0; i < d1; i++) {
1149 #pragma _CRI concurrent
1150 for (j = 0; j < d2; j++) {
1152 if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1153 if (cell_samples[hr][pos] == 0)
1154 means[hr][pos] = val[pos];
1156 means[hr][pos] += val[pos];
1157 ++cell_samples[hr][pos];
1166 void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
1170 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1171 for (i = 0; i < d1; i++) {
1172 #pragma _CRI concurrent
1173 for (j = 0; j < d2; j++) {
1175 if (cell_samples[hr][pos] != 0) {
1176 means[hr][pos] /= (double)cell_samples[hr][pos];
1185 void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
1186 char FillFlag, float FillValue, double **means, double **stddevs,
1187 size_t **cell_samples)
1192 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1194 for (i = 0; i < d1; i++) {
1195 #pragma _CRI concurrent
1196 for (j = 0; j < d2; j++) {
1198 if (!FillFlag || (FillFlag && val[pos] != FillValue
1199 && means[hr][pos] != FillValue)) {
1200 delta = means[hr][pos] - (double)val[pos];
1201 if (cell_samples[hr][pos] == 0)
1202 stddevs[hr][pos] = delta * delta;
1204 stddevs[hr][pos] += delta * delta;
1205 ++cell_samples[hr][pos];
1213 void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
1214 char FillFlag, float FillValue, double **means, double **stddevs,
1215 size_t **cell_samples)
1220 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1222 for (i = 0; i < d1; i++) {
1223 #pragma _CRI concurrent
1224 for (j = 0; j < d2; j++) {
1226 if (!FillFlag || (FillFlag && val[pos] != FillValue
1227 && means[hr][pos] != FillValue)) {
1228 delta = means[hr][pos] - val[pos];
1229 if (cell_samples[hr][pos] == 0)
1230 stddevs[hr][pos] = delta * delta;
1232 stddevs[hr][pos] += delta * delta;
1233 ++cell_samples[hr][pos];
1242 void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
1246 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1247 for (i = 0; i < d1; i++) {
1248 #pragma _CRI concurrent
1249 for (j = 0; j < d2; j++) {
1251 if (cell_samples[hr][pos] != 0) {
1252 stddevs[hr][pos] /= (double)cell_samples[hr][pos];
1253 stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
1263 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)
1267 size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1268 char name[NC_MAX_NAME+1];
1271 /* Allocate space for timeslice */
1272 val = alloc_var(nctype, (d1 * d2));
1274 for (i = 0; i < nifnames; i++) {
1276 printf("\tOpening %s", input_fnames[i]);
1277 if (ndims > 3) printf(" to retrieve level %d", level);
1278 if (stat_type == MEAN_TYPE)
1279 printf(" for computing means\n");
1281 printf(" for computing stddevs\n");
1283 /* Open input file */
1284 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1286 * Inquire about number of dimensions, variables, global
1287 * attributes and the ID of the unlimited dimension
1289 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1290 &input_ngatts, &input_unlimdimid));
1291 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1292 if (strcmp(name, time_name)) {
1293 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1296 /* Get seconds of day variable */
1297 wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
1298 if (!(mcsec = (int *)malloc(len * sizeof(int)))) {
1299 perror("read_and_compute:mcsec");
1302 wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));
1303 /* Get variable ID for requested variable */
1304 wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
1305 /* Retrieve time slice of desired variable */
1306 for (ts = 0; ts < len; ts++) {
1309 start[(ndims-2)] = 0;
1310 count[(ndims-2)] = d1;
1311 start[(ndims-1)] = 0;
1312 count[(ndims-1)] = d2;
1319 wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));
1320 if (stat_type == MEAN_TYPE)
1321 add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1323 add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1326 wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));
1327 if (stat_type == MEAN_TYPE)
1328 add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1330 add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1333 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1338 /* Free mcsec vector */
1340 /* Close input file */
1341 wrap_nc(nc_close(input_ncid));
1343 /* Divide sums by number of samples */
1344 if (stat_type == MEAN_TYPE)
1345 divide_means(d1, d2, means, cell_samples);
1347 divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
1349 /* Free working variable array */
1355 float *double_to_float(size_t len, double *dvar)
1360 if (!(fvar = (float *)malloc(len * sizeof(float)))) {
1361 perror("double_to_float:fvar");
1365 for (i = 0; i < len; i++)
1366 fvar[i] = (float)dvar[i];
1371 void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
1372 int level, int ndims, double **var_hours)
1374 /* Output dimensions should be one larger than input dimensions */
1376 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1379 if (nctype == NC_FLOAT) {
1380 if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
1381 perror("write_var_hours:fvar");
1386 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1388 count[0] = 1; /* time */
1390 count[1] = 1; /* hour */
1391 start[(ndims-2)] = 0;
1392 count[(ndims-2)] = d1;
1393 start[(ndims-1)] = 0;
1394 count[(ndims-1)] = d2;
1401 for (i = 0; i < (d1 * d2); i++)
1402 fvar[i] = (float)var_hours[hr][i];
1403 wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
1406 wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
1409 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1414 if (nctype == NC_FLOAT)
1420 void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
1423 size_t d1, d2, **cell_samples;
1424 double **means, **stddevs;
1425 struct var *in_vnode, *out_vnode;
1427 if (input_var_head) {
1428 in_vnode = input_var_head;
1429 /* March through non-metadata variables to compute stats */
1430 for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
1431 if (!in_vnode->metadata) {
1432 /* Make sure time is really the first dimension */
1433 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
1434 fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
1437 /* Figure out input dimensions */
1438 if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
1439 fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
1442 /* Assume 2D output; find dimensions */
1443 d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
1444 d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
1445 if (in_vnode->ndims == 3)
1448 nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
1449 /* Allocate working space for means and stddevs */
1450 alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
1451 init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
1452 printf("Computing statistics for %s\n",
1454 /* Compute means and stddevs, then write them */
1455 for (k = 0; k < nlevels; k++) {
1456 /* Read and compute means */
1457 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);
1458 /* Find corresponding output variable on the mean output file */
1459 out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
1460 /* Write out the means for this variable */
1461 write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
1462 /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
1463 reset_cell_samples(d1, d2, cell_samples);
1464 /* Read and compute stddevs using means */
1465 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);
1466 /* Find corresponding output variable on the stddev output file */
1467 out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
1468 /* Write out stddevs for this variable */
1469 write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
1472 /* Free working space for means and stddevs */
1473 free_means_stddevs(means, stddevs, cell_samples);
1475 in_vnode = in_vnode->next;
1481 void usage(char *program)
1483 fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
1484 fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
1488 int main(int argc, char **argv)
1490 int i, ic, nifnames;
1491 size_t len, pos, nsamples;
1492 char *mfname, *sfname, **ifnames, *flist;
1495 mfname = sfname = flist = NULL;
1499 printf("Number of metadata variables (nmvars) = %d\n", nmvars);
1503 /* check command-line flags and switches */
1504 while ((ic = getopt(argc, argv, "m:s:")) != -1) {
1507 if (!(mfname = strdup(optarg))) {
1513 if (!(sfname = strdup(optarg))) {
1522 fprintf(stderr, "Output file name for writing means required.\n");
1526 fprintf(stderr, "Output file name for writing standard deviations required.\n");
1530 if (optind < argc) {
1531 if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
1535 for (i = optind; i < argc; i++) {
1536 if (!(ifnames[nifnames++] = strdup(argv[i]))) {
1537 perror("ifnames[nifnames++]");
1543 fprintf(stderr, "At least one input file name is required.\n");
1549 * Build list of input files to be included in the global history
1550 * attribute on the two outputs files.
1552 for (i = len = 0; i < nifnames; i++)
1553 len += strlen(ifnames[i]);
1554 len += nifnames + 1; /* space in front of every name + null terminator */
1555 if (!(flist = (char *)malloc(len * sizeof(char)))) {
1559 for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
1560 sprintf(&flist[pos], " %s", ifnames[i]);
1562 printf("flist=%s\n", flist);
1566 /* Open every file to count up the number of time samples. */
1567 nsamples = count_samples(nifnames, ifnames);
1569 printf("Number of samples across all files = %ld\n", (long)nsamples);
1574 * Open the *last* input file and the two output files (for mean and
1575 * standard deviation). Define dimensions, variables, and attributes
1576 * for the two output files based on the *last* input files. The
1577 * last file is used because some metadata variables will contain
1578 * only values from the last time slice from the period over which
1579 * the statistics are computed.
1581 open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
1584 compute_stats(nifnames, ifnames, nsamples);
1586 /* Close the two output files */
1587 wrap_nc(nc_close(mean_ncid));
1588 wrap_nc(nc_close(stddev_ncid));