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]);
1262 void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp)
1267 /* Allocate space for all timeslices */
1268 if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) {
1269 perror("alloc_all_samples:*valsp");
1273 for (i = 0; i < nsamples; i++) {
1274 vals[i] = alloc_var(nctype, (d1 * d2));
1276 if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) {
1277 perror("alloc_all_samples:*mcsecp");
1284 void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec)
1287 float **fvals = NULL;
1288 double **dvals = NULL;
1292 fvals = (float **)vals;
1295 dvals = (double **)vals;
1298 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1302 for (i = 0; i < nsamples; i++) {
1303 for (j = 0; j < (d1 * d2); j++) {
1306 fvals[i][j] = (float)0.0;
1309 dvals[i][j] = (double)0.0;
1312 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1322 void free_all_samples(size_t nsamples, void **vals, int *mcsec)
1326 for (i = 0; i < nsamples; i++)
1327 if (vals[i]) free(vals[i]);
1335 void read_all_samples(int nifnames, char **input_fnames, size_t d1, size_t d2, size_t nsamples, char *var_name, nc_type nctype, int level, int ndims, void **vals, int *mcsec)
1339 size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1340 char name[NC_MAX_NAME+1];
1342 for (i = 0; i < nifnames; i++) {
1344 printf("\tOpening %s", input_fnames[i]);
1345 if (ndims > 3) printf(" to retrieve level %d\n", level);
1347 /* Open input file */
1348 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1350 * Inquire about number of dimensions, variables, global
1351 * attributes and the ID of the unlimited dimension
1353 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1354 &input_ngatts, &input_unlimdimid));
1355 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1356 if (strcmp(name, time_name)) {
1357 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1360 /* Make sure we don't run off the end of allocated space */
1361 if ((nslice+len) > nsamples) {
1362 fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);
1365 /* Get seconds of day variable */
1366 wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
1367 wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));
1368 /* Get variable ID for requested variable */
1369 wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
1370 /* Retrieve time slice of desired variable */
1371 for (ts = 0; ts < len; ts++) {
1374 start[(ndims-2)] = 0;
1375 count[(ndims-2)] = d1;
1376 start[(ndims-1)] = 0;
1377 count[(ndims-1)] = d2;
1384 wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));
1387 wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));
1390 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1396 /* Close input file */
1397 wrap_nc(nc_close(input_ncid));
1403 void compute(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples)
1407 for (i = 0; i < nsamples; i++) {
1410 if (stat_type == MEAN_TYPE)
1411 add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
1413 add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1416 if (stat_type == MEAN_TYPE)
1417 add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
1419 add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1422 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1427 /* Divide sums by number of samples */
1428 if (stat_type == MEAN_TYPE)
1429 divide_means(d1, d2, means, cell_samples);
1431 divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
1436 float *double_to_float(size_t len, double *dvar)
1441 if (!(fvar = (float *)malloc(len * sizeof(float)))) {
1442 perror("double_to_float:fvar");
1446 for (i = 0; i < len; i++)
1447 fvar[i] = (float)dvar[i];
1452 void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
1453 int level, int ndims, double **var_hours)
1455 /* Output dimensions should be one larger than input dimensions */
1457 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1460 if (nctype == NC_FLOAT) {
1461 if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
1462 perror("write_var_hours:fvar");
1467 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1469 count[0] = 1; /* time */
1471 count[1] = 1; /* hour */
1472 start[(ndims-2)] = 0;
1473 count[(ndims-2)] = d1;
1474 start[(ndims-1)] = 0;
1475 count[(ndims-1)] = d2;
1482 for (i = 0; i < (d1 * d2); i++)
1483 fvar[i] = (float)var_hours[hr][i];
1484 wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
1487 wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
1490 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1495 if (nctype == NC_FLOAT)
1501 void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
1503 int j, k, nlevels, *mcsec;
1504 size_t d1, d2, **cell_samples;
1505 double **means, **stddevs;
1506 struct var *in_vnode, *out_vnode;
1509 if (input_var_head) {
1510 in_vnode = input_var_head;
1511 /* March through non-metadata variables to compute stats */
1512 for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
1513 if (!in_vnode->metadata) {
1514 /* Make sure time is really the first dimension */
1515 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
1516 fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
1519 /* Figure out input dimensions */
1520 if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
1521 fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
1524 /* Assume 2D output; find dimensions */
1525 d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
1526 d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
1527 if (in_vnode->ndims == 3)
1530 nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
1531 /* Allocate working space for means and stddevs */
1532 alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
1533 init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
1534 /* Allocate and initize space for entire field across all files */
1535 alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);
1536 init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);
1537 printf("Computing statistics for %s\n",
1539 /* Compute means and stddevs, then write them */
1540 for (k = 0; k < nlevels; k++) {
1541 /* Read the entire fields from all files */
1542 read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);
1544 compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
1545 /* Find corresponding output variable on the mean output file */
1546 out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
1547 /* Write out the means for this variable */
1548 write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
1549 /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
1550 reset_cell_samples(d1, d2, cell_samples);
1551 /* Compute stddevs using means */
1552 compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
1553 /* Find corresponding output variable on the stddev output file */
1554 out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
1555 /* Write out stddevs for this variable */
1556 write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
1559 /* Free all samples */
1560 free_all_samples(nsamples, vals, mcsec);
1561 /* Free working space for means and stddevs */
1562 free_means_stddevs(means, stddevs, cell_samples);
1564 in_vnode = in_vnode->next;
1570 void usage(char *program)
1572 fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
1573 fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
1577 int main(int argc, char **argv)
1579 int i, ic, nifnames;
1580 size_t len, pos, nsamples;
1581 char *mfname, *sfname, **ifnames, *flist;
1584 mfname = sfname = flist = NULL;
1588 printf("Number of metadata variables (nmvars) = %d\n", nmvars);
1592 /* check command-line flags and switches */
1593 while ((ic = getopt(argc, argv, "m:s:")) != -1) {
1596 if (!(mfname = strdup(optarg))) {
1602 if (!(sfname = strdup(optarg))) {
1611 fprintf(stderr, "Output file name for writing means required.\n");
1615 fprintf(stderr, "Output file name for writing standard deviations required.\n");
1619 if (optind < argc) {
1620 if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
1624 for (i = optind; i < argc; i++) {
1625 if (!(ifnames[nifnames++] = strdup(argv[i]))) {
1626 perror("ifnames[nifnames++]");
1632 fprintf(stderr, "At least one input file name is required.\n");
1638 * Build list of input files to be included in the global history
1639 * attribute on the two outputs files.
1641 for (i = len = 0; i < nifnames; i++)
1642 len += strlen(ifnames[i]);
1643 len += nifnames + 1; /* space in front of every name + null terminator */
1644 if (!(flist = (char *)malloc(len * sizeof(char)))) {
1648 for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
1649 sprintf(&flist[pos], " %s", ifnames[i]);
1651 printf("flist=%s\n", flist);
1655 /* Open every file to count up the number of time samples. */
1656 nsamples = count_samples(nifnames, ifnames);
1658 printf("Number of samples across all files = %ld\n", (long)nsamples);
1663 * Open the *last* input file and the two output files (for mean and
1664 * standard deviation). Define dimensions, variables, and attributes
1665 * for the two output files based on the *last* input files. The
1666 * last file is used because some metadata variables will contain
1667 * only values from the last time slice from the period over which
1668 * the statistics are computed.
1670 open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
1673 compute_stats(nifnames, ifnames, nsamples);
1675 /* Close the two output files */
1676 wrap_nc(nc_close(mean_ncid));
1677 wrap_nc(nc_close(stddev_ncid));