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
18 #define HOURS_PER_DAY 24
19 #define MIN_PER_HOUR 60
20 #define SEC_PER_MIN 60
21 #define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
23 static char *metadata_vars[] = {
29 "hour", /* new metadata variable to be added to output files */
59 struct text_att *next;
60 struct text_att *prev;
85 static char *time_name = "time";
86 static char *time_bounds_attname = "bounds";
87 static char *time_bounds_name = "time_bounds";
88 static char *climatology_bounds_attname = "climatology";
89 static char *climatology_bounds_name = "climatology_bounds";
90 static char *mcdate_name = "mcdate";
91 static char *mcsec_name = "mcsec";
92 static char *history_name = "history";
93 static char *nsamples_name = "nsamples";
94 static char *cell_method_name = "cell_method";
95 static char *cell_methods_prefix[] = { "time: mean within days time: mean over days", "time: mean within days times: standard_deviation over days" };
97 static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
98 static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
99 static struct text_att *input_gatt_head = NULL;
100 static struct dim *input_dim_head = NULL, **input_dim_idx;
101 static struct var *input_var_head = NULL;
102 /* translation stuff */
103 static int *idid2mdid, *idid2sdid; /* dimension IDs */
105 static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
106 static struct dim *mean_dim_head = NULL;
107 static struct var *mean_var_head = NULL;
108 static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
109 static struct dim *stddev_dim_head = NULL;
110 static struct var *stddev_var_head = NULL;
112 static double *times;
113 static double *tbounds;
117 char is_metadata(char *name)
121 for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
130 struct dim *dimlist_find_by_name(struct dim *head, char *name)
133 struct dim *p = NULL;
137 for (i = 0 ; i == 0 || node != head; i++) {
138 if (!strcmp(node->name, name))
149 struct var *varlist_find_by_name(struct var *head, char *name)
156 for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
158 if (i && node == head)
167 void gattlist_add_node(struct text_att **headp, char *name, char *value)
169 struct text_att *head, *node;
173 if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
174 perror("gattlist_add_node");
178 if (!(node->name = strdup(name))) {
179 perror("gattlist_add_node");
182 if (!(node->value = strdup(value))) {
183 perror("gattlist_add_node");
195 node->prev = head->prev;
196 /* set this after setting node->prev from here */
198 /* set this after having set node->prev */
199 node->prev->next = node;
205 struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
207 struct dim *head, *node;
211 if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
212 perror("dimlist_add_node");
217 if (!(node->name = strdup(name))) {
218 perror("dimlist_add_node");
231 node->prev = head->prev;
232 /* set this after setting node->prev from here */
234 /* set this after having set node->prev */
235 node->prev->next = node;
241 void varlist_add_node(struct var **headp, int ncvarid, char *name,
242 nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
246 struct var *head, *node;
250 if (!(node = (struct var *)malloc(sizeof(struct var)))) {
251 perror("varlist_add_node");
255 node->ncvarid = ncvarid;
256 if (!(node->name = strdup(name))) {
257 perror("varlist_add_node");
260 node->nctype = nctype;
262 if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
263 perror("varlist_add_node");
266 for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
268 node->metadata = is_metadata(name);
269 node->FillFlag = FillFlag;
270 node->FillValue = FillValue;
280 node->prev = head->prev;
281 /* set this after setting node->prev from here */
283 /* set this after having set node->prev */
284 node->prev->next = node;
290 void varlist_print(struct var *headp)
295 printf("Beginning of Variable List\n");
299 for (j = 0; j == 0 || node != headp; j++) {
300 printf("Variable %d (ptr=%ld):\n", j, (long)node);
301 printf("\tncvarid = %d\n", node->ncvarid);
302 printf("\tname = %s\n", node->name);
303 printf("\tnctype = %d\n", node->nctype);
304 printf("\tndims = %d\n", node->ndims);
305 printf("\tdimids =");
306 for (i = 0; i < node->ndims; i++)
307 printf(" %d", node->dimids[i]);
309 printf("\tnatts = %d\n", node->natts);
310 printf("\tmetadata = %d\n", (int)node->metadata);
311 printf("\tnext = %ld\n", (long)node->next);
312 printf("\tprev = %ld\n", (long)node->prev);
317 printf("\tList undefined\n");
320 printf("End of Variable List\n");
325 void wrap_nc(int status)
327 if (status != NC_NOERR) {
328 fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
335 void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
338 char name[NC_MAX_NAME+1];
340 struct dim **in_dim_idx;
342 if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
343 perror("*in_dim_idxp");
346 in_dim_idx = *in_dim_idxp;
348 for (i = 0; i < ndims; i++) {
349 wrap_nc(nc_inq_dim(ncid, i, name, &len));
350 in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
355 void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
358 char name[NC_MAX_NAME+1], *value;
362 for (i = 0; i < ngatts; i++) {
363 wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
364 wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
365 if (xtype != NC_CHAR) {
366 fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
369 if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
370 perror("get_gatts: value");
373 wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
374 value[(len+1-1)] = '\0';
375 gattlist_add_node(gatt_headp, name, value);
376 free(value); /* because gattlist_add_node() duplicates it */
382 void get_vars(int ncid, int nvars, struct var **var_headp)
384 int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
387 char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
388 nc_type xtype, atype;
390 for (i = 0; i < nvars; i++) {
391 wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
393 /* Check for _FillValue */
396 if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
398 if (atype == NC_FLOAT && len) {
399 wrap_nc(nc_get_att_float(ncid, i,
400 fill_att_name, &FillValue));
404 varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
405 natts, FillFlag, FillValue);
411 int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
412 size_t nsamples, int **in2outp, int out_ncid,
413 struct dim **out_dim_headp, int *out_unlimdimidp)
416 * Define dimensions on new files and build translation tables between
419 int j, dimid, ndims, *in2out;
425 if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
426 perror("put_dimensions:in2outp");
433 for (j = 0; j == 0 || dnode != in_dim_head; j++) {
434 if (dnode->dimid == in_unlimdimid)
438 wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
439 dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
441 in2out[dnode->dimid] = dimid;
442 if (dnode->dimid == in_unlimdimid)
443 *out_unlimdimidp = dimid;
444 /* Just after the "time" dimension, add the new
445 * "nsamples" dimension. */
446 if (!strcmp(dnode->name, time_name)) {
447 wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
448 dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
455 fprintf(stderr, "WARNING: No dimensions defined!\n");
461 int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
464 * Write out global attributes matching those of the input file.
465 * Change history attribute to the string passed in as an argument.
467 int j, hflag = 0, ngatts;
469 struct text_att *anode;
474 anode = in_gatt_head;
475 for (j = 0; j == 0 || anode != in_gatt_head; j++) {
476 if (!strcmp(anode->name, history_name)) {
481 value = anode->value;
482 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
486 /* If no history attribute on input, add one to the output */
488 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
493 fprintf(stderr, "WARNING: No global attributes defined!\n");
499 int translate_dimids(int in_ndims, int *in_dimids, int *in2out, int *out_dimids)
502 * Translate between input dimension IDs and those for the output file.
506 for (i = ndims = 0; i < in_ndims; i++, ndims++)
507 out_dimids[ndims] = in2out[in_dimids[i]];
512 int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
513 int in_varid, int out_ncid, int out_varid, char **cell_methods)
516 * Copy the attributes of the input variable to those of the output
517 * variable. If the variable is not a metadata variable, ensure that
518 * the cell_method attribute is set correctly, depending on the output
523 char name[NC_MAX_NAME+1], cmflag = 0;
525 for (i = natts = 0; i < in_natts; i++, natts++) {
526 wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
527 if (!strcmp(name, cell_method_name) && !metadata) {
528 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
531 else if (!strcmp(name, time_bounds_attname))
532 /* Change time_bounds to climatology_bounds */
533 wrap_nc(nc_put_att_text(out_ncid, out_varid, climatology_bounds_attname, strlen(climatology_bounds_name), climatology_bounds_name));
535 wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
538 * If no cell_method attribute was specified for a non-metadata
539 * variable on the input file, add the appropriate cell_method anyway
540 * on the output file.
542 if (!cmflag && !metadata) {
543 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
550 int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
551 int *in2out, int out_ncid, struct var **out_var_headp,
552 char **cell_methods, int stat_type)
555 * Define variables on output file and copy attributes from input file.
556 * Return number of variables defined.
558 int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
566 * March through input variables creating (mostly) the same
567 * variables on the output file.
569 for (j = 0; j == 0 || vnode != in_var_head; j++) {
570 ndims = translate_dimids(vnode->ndims, vnode->dimids, in2out, dimids);
571 /* Create all normal data variables, but only the four
572 * time-based metadata variables: time,
573 * climatology_bounds, mcdate, and mcsec. */
574 if (!strcmp(vnode->name, time_name)) {
575 /* Magically promote "time" variable to
576 * double precision on output file */
577 wrap_nc(nc_def_var(out_ncid, vnode->name, NC_DOUBLE, ndims, dimids, &varid));
578 natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);
579 varlist_add_node(out_var_headp, varid, vnode->name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
581 /* Create "climatology_bounds" right after
582 * time and instead of "time_bounds", with no
584 wrap_nc(nc_def_var(out_ncid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, &varid));
586 varlist_add_node(out_var_headp, varid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
589 else if (!vnode->metadata || !strcmp(vnode->name, mcdate_name) || !strcmp(vnode->name, mcsec_name) || (vnode->metadata && strcmp((in_dim_idx[vnode->dimids[0]])->name, time_name))) {
590 /* If it is not a metadata variable OR it is
591 * mcdate OR it is mcsec OR (it is a metadata
592 * variable that does not have time as its
593 * first dimension), then create it */
594 wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
595 natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);
596 varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
603 fprintf(stderr, "WARNING: No variables defined!\n");
609 void *alloc_var(nc_type nctype, size_t len)
615 if (!(val = (float *)malloc((len) * sizeof(float)))) {
616 perror("alloc_var: val");
621 if (!(val = (int *)malloc((len) * sizeof(int)))) {
622 perror("alloc_var: val");
627 if (!(val = (double *)malloc((len) * sizeof(double)))) {
628 perror("alloc_var: val");
633 if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
634 perror("alloc_var: val");
639 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
646 void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
649 size_t len = (size_t)1;
652 /* Figure out the total size */
653 for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
655 val = alloc_var(nctype, len);
659 wrap_nc(nc_get_var_float(ncid, varid, val));
662 wrap_nc(nc_get_var_int(ncid, varid, val));
665 wrap_nc(nc_get_var_double(ncid, varid, val));
668 wrap_nc(nc_get_var_text(ncid, varid, val));
671 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
678 void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
681 size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
684 /* Make sure time is really the first dimension */
685 if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
686 fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
690 * Figure out the total size for the slice (start at second dimension)
691 * and build start/count vectors.
695 for (i = 1; i < ndims; i++) {
697 count[i] = (dim_idx[dimids[i]])->len;
698 len *= (dim_idx[dimids[i]])->len;
701 val = alloc_var(nctype, len);
705 wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
708 wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
711 wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
714 wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
717 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
724 void write_var(int ncid, int varid, nc_type nctype, void *val)
728 wrap_nc(nc_put_var_float(ncid, varid, val));
731 wrap_nc(nc_put_var_int(ncid, varid, val));
734 wrap_nc(nc_put_var_double(ncid, varid, val));
737 wrap_nc(nc_put_var_text(ncid, varid, val));
740 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
747 void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
750 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
752 /* Make sure time is really the first dimension */
753 if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
754 fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
758 /* Build start/count vectors */
761 for (i = 1; i < ndims; i++) {
763 count[i] = (dim_idx[dimids[i]])->len;
768 wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
771 wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
774 wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
777 wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
780 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
787 void write_time_var(int ncid, int varid, nc_type nctype, void *time_var)
789 size_t start[1], count[1];
796 wrap_nc(nc_put_vara_float(ncid, varid, start, count, time_var));
799 wrap_nc(nc_put_vara_int(ncid, varid, start, count, time_var));
802 wrap_nc(nc_put_vara_double(ncid, varid, start, count, time_var));
805 wrap_nc(nc_put_vara_text(ncid, varid, start, count, time_var));
808 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
815 void write_tbounds(int ncid, int varid, double *tbounds)
817 size_t start[2], count[2];
819 start[0] = start[1] = 0;
820 count[0] = 24; count[1] = 2;
822 wrap_nc(nc_put_vara_double(ncid, varid, start, count, tbounds));
827 void copy_metadata(int in_ncid, struct var *in_var_head,
828 struct dim **in_dim_idx, int out_ncid, struct var *out_var_head,
829 double *times, double *tbounds, int *mcdate, int *mcsec)
832 struct var *in_vnode, *out_vnode;
835 if (!in_var_head) return;
837 in_vnode = in_var_head;
839 * March through input variables to stuff metadata values into
840 * the new output files. NOTE: The only time-based metadata
841 * variables, explicitly handled, are: time, climatology_bounds,
842 * mcdate, and mcsec. These are passed in.
844 for (i = 0; i == 0 || in_vnode != in_var_head; i++) {
845 if (in_vnode->metadata) {
846 /* Find the corresponding output variable */
847 if (!strcmp(in_vnode->name, time_bounds_name))
848 out_vnode = varlist_find_by_name(out_var_head, climatology_bounds_name);
850 out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
852 /* Read values from input and write to output */
853 if (strcmp((in_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
854 /* "time" is not the first dimension,
855 * so just copy the variable values */
857 printf("Copying metadata variable %s\n",
860 val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
861 write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
865 /* "time" is the first dimension, so
866 * explicitly handle the "time",
867 * "climatology_bounds", "mcdate", and
868 * "mcsec" while dropping all others */
869 if (!strcmp(out_vnode->name, time_name)) {
870 /* Force time stamps to be
871 * those computed previously,
872 * and remember these are now
873 * double precision. */
875 printf("Writing metadata variable %s\n", out_vnode->name);
877 write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, times);
879 if (!strcmp(out_vnode->name, climatology_bounds_name)) {
880 /* time_bounds are now
881 * climatology_bounds and have
882 * pre-computed values */
884 printf("Writing metadata variable %s\n", out_vnode->name);
886 write_tbounds(out_ncid, out_vnode->ncvarid, tbounds);
888 else if (!strcmp(out_vnode->name, mcdate_name)) {
889 /* mcdate is pre-computed */
891 printf("Writing metadata variable %s\n", out_vnode->name);
893 write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcdate);
895 else if (!strcmp(out_vnode->name, mcsec_name)) {
896 /* mcsec is pre-computed */
898 printf("Writing metadata variable %s\n", out_vnode->name);
900 write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcsec);
906 printf("Dropping metadata variable %s\n", in_vnode->name);
912 printf("Skipping non-metadata variable %s\n",
919 in_vnode = in_vnode->next;
925 void get_time_bounds(char *first_fname, char *last_fname, double *times, double *tbounds, int *mcdate, int *mcsec)
927 int i, time_dimid, time_bounds_varid, mcdate_varid, mcsec_varid;
928 size_t time_len, time_bounds_start[2], time_bounds_count[2],
929 mc_start[1], mc_count[1];
930 double bnd1[24], bnd2[24], day1, day2;
932 /* Open first input file */
933 wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
934 /* Get dimension ID of the time dimension */
935 wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
936 /* Get length of time dimension */
937 wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
938 /* Get variable ID of time_bounds variable */
939 wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
940 /* Set start/count for reading the left values of time_bounds from the
941 * first 24 timeslices of the first file */
942 time_bounds_start[0] = time_bounds_start[1] = 0;
943 time_bounds_count[0] = 24; time_bounds_count[1] = 1;
944 /* Read the values */
945 wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,
946 time_bounds_start, time_bounds_count, bnd1));
947 /* If the first and last file are not the same, close the first one and
948 * open the second one */
949 if (strcmp(first_fname, last_fname)) {
950 /* Close the first input file */
951 wrap_nc(nc_close(input_ncid));
952 /* Open the last input file */
953 wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
954 /* Get dimension ID of the time dimension */
955 wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
956 /* Get length of time dimension */
957 wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
958 /* Get variable ID of time_bounds variable */
959 wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
960 &time_bounds_varid));
962 /* Set start/count for reading the right values of time_bounds from the
963 * last 24 timeslices of the last file */
964 time_bounds_start[0] = time_len - 24; time_bounds_start[1] = 1;
965 time_bounds_count[0] = 24; time_bounds_count[1] = 1;
967 wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,
968 time_bounds_start, time_bounds_count, bnd2));
969 /* Read the last 24 mcdate and mcsec values */
970 wrap_nc(nc_inq_varid(input_ncid, mcdate_name, &mcdate_varid));
971 wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &mcsec_varid));
972 mc_start[0] = time_len - 24;
974 wrap_nc(nc_get_vara_int(input_ncid, mcdate_varid, mc_start, mc_count,
976 wrap_nc(nc_get_vara_int(input_ncid, mcsec_varid, mc_start, mc_count,
979 /* Close the last input file */
980 wrap_nc(nc_close(input_ncid));
982 /* Store time bounds and compute new time stamps as midpoints */
983 for (i = 0; i < 24; i++) {
984 tbounds[(i*2)] = bnd1[i];
985 tbounds[(i*2+1)] = bnd2[i];
986 /* Pick one (integer) day near the midpoint of the bounds and
987 * add in the fraction from the right bounds, so for January
988 * the first hour-of-day value would be 15.0 (plus some year
990 day1 = floor(bnd1[i]);
991 day2 = floor(bnd2[i]);
992 times[i] = floor((day1 + day2) / 2.0) + (bnd2[i]-day2);
998 char **make_cell_methods(size_t nsamples)
1002 char **cell_methods;
1004 if (!(cell_methods = (char **)malloc((NUM_TYPES * sizeof(char *))))) {
1005 perror("make_cell_methods:cell_methods");
1008 for (i = 0; i < NUM_TYPES; i++) {
1009 len = strlen(cell_methods_prefix[i]) + 48;
1010 if (!(cell_methods[i] = (char *)malloc((len * sizeof(char))))) {
1011 perror("make_cell_methods:cell_methods[i]");
1014 sprintf(cell_methods[i], "%s (comment: %ld samples)",
1015 cell_methods_prefix[i], (long)nsamples);
1017 return cell_methods;
1020 void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
1023 char *mean_history_gatt, *stddev_history_gatt,
1024 *mean_prefix="Statistical means from history files:",
1025 *stddev_prefix="Statistical standard deviations from history files:",
1029 * Construct strings for history global attributes for the two output
1032 if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
1033 perror("open_inout:mean_history_gatt");
1036 if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
1037 perror("open_inout:stddev_history_gatt");
1040 sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
1041 sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
1043 /* Allocate space for time values for each hour of the day */
1044 if (!(times = (double *)malloc((24 * sizeof(double))))) {
1045 perror("open_inout:times");
1048 /* Allocate space for two time bounds values for each hour of the day */
1049 if (!(tbounds = (double *)malloc((24 * 2 * sizeof(double))))) {
1050 perror("open_inout:tbounds");
1053 /* Allocate space for current date and seconds for each hour of day */
1054 if (!(mcdate = (int *)malloc((24 * sizeof(int))))) {
1055 perror("open_inout:mcdate");
1058 if (!(mcsec = (int *)malloc((24 * sizeof(int))))) {
1059 perror("open_inout:mcsec");
1065 * Explicitly handle "time", "time_bounds" -> "climatology_bounds",
1066 * mcdate, and mcsec metadata variables by constructing their values
1067 * before doing anything else.
1069 get_time_bounds(first_fname, last_fname, times, tbounds, mcdate, mcsec);
1071 for (i = 0; i < 24; i++)
1072 fprintf(stderr, "Hour of day=%d, time=%lf, tbounds=%lf,%lf mcdate=%d mcsec=%d\n", i, times[i], tbounds[(i*2)], tbounds[(i*2+1)], mcdate[i], mcsec[i]);
1075 /* Open last input file */
1076 wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
1077 /* Inquire about number of dimensions, variables, global attributes
1078 * and the ID of the unlimited dimension
1080 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
1081 &input_unlimdimid));
1082 /* Grab dimension IDs and lengths from input file */
1083 get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
1084 /* Grab desired global attributes from input file */
1085 get_gatts(input_ncid, input_ngatts, &input_gatt_head);
1086 /* Grab variable IDs and attributes from input file */
1087 get_vars(input_ncid, input_nvars, &input_var_head);
1089 varlist_print(input_var_head);
1091 /* Create netCDF files for monthly means and standard deviations */
1092 /* Create new netCDF files */
1093 wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
1094 wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
1095 /* Define dimensions */
1096 mean_ndims = put_dimensions(input_dim_head, input_ndims,
1097 input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
1098 &mean_dim_head, &mean_unlimdimid);
1099 stddev_ndims = put_dimensions(input_dim_head, input_ndims,
1100 input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
1101 &stddev_dim_head, &stddev_unlimdimid);
1102 /* Define cell_methods[], including number of samples */
1103 cell_methods = make_cell_methods(nsamples);
1104 /* Define variables and their attributes */
1105 mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
1106 idid2mdid, mean_ncid, &mean_var_head, cell_methods,
1108 stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
1109 idid2sdid, stddev_ncid, &stddev_var_head, cell_methods,
1112 printf("Mean variables:\n");
1113 varlist_print(mean_var_head);
1114 printf("Stddev variables:\n");
1115 varlist_print(stddev_var_head);
1117 /* Store global attributes */
1118 mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
1119 stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
1120 stddev_history_gatt);
1121 /* End define mode */
1122 wrap_nc(nc_enddef(mean_ncid));
1123 wrap_nc(nc_enddef(stddev_ncid));
1124 /* Write out metdata variables that do not depend on "time" */
1125 copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
1126 mean_var_head, times, tbounds, mcdate, mcsec);
1127 copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
1128 stddev_var_head, times, tbounds, mcdate, mcsec);
1130 wrap_nc(nc_close(input_ncid));
1135 size_t count_samples(int nifnames, char **input_fnames)
1138 char name[NC_MAX_NAME+1];
1139 size_t len, cnt = 0;
1141 for (i = 0; i < nifnames; i++) {
1142 /* Open input file */
1143 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1145 * Inquire about number of dimensions, variables, global
1146 * attributes and the ID of the unlimited dimension
1148 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1149 &input_ngatts, &input_unlimdimid));
1150 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1151 if (strcmp(name, time_name)) {
1152 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1156 printf("%ld samples in %s\n", (long)len, input_fnames[i]);
1158 wrap_nc(nc_close(input_ncid));
1164 void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
1167 * Allocate space for arrays of means and standard deviations by
1171 size_t **cell_samples;
1172 double **means, **stddevs;
1174 if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1175 perror("alloc_means_stddevs:*meansp");
1178 if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1179 perror("alloc_means_stddevs:*stddevsp");
1182 if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
1183 perror("alloc_means_stddevs:*cell_samplesp");
1188 stddevs = *stddevsp;
1189 cell_samples = *cell_samplesp;
1191 for (i = 0; i < HOURS_PER_DAY; i++) {
1192 if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1193 perror("alloc_means_stddevs:means[i]");
1196 if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1197 perror("alloc_means_stddevs:stddevs[i]");
1200 if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
1201 perror("alloc_means_stddevs:cell_samples[i]");
1209 void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
1213 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1214 for (i = 0; i < d1; i++) {
1215 #pragma _CRI concurrent
1216 for (j = 0; j < d2; j++) {
1218 means[hr][pos] = FillValue;
1219 stddevs[hr][pos] = FillValue;
1220 cell_samples[hr][pos] = 0;
1228 void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
1232 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1233 for (i = 0; i < d1; i++) {
1234 #pragma _CRI concurrent
1235 for (j = 0; j < d2; j++) {
1237 cell_samples[hr][pos] = 0;
1245 void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
1248 * Free space from arrays of means and standard deviations by
1253 for (i = 0; i < HOURS_PER_DAY; i++) {
1256 free(cell_samples[i]);
1266 void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
1267 char FillFlag, float FillValue, double **means, size_t **cell_samples)
1271 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1273 for (i = 0; i < d1; i++) {
1274 #pragma _CRI concurrent
1275 for (j = 0; j < d2; j++) {
1277 if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1278 if (cell_samples[hr][pos] == 0)
1279 means[hr][pos] = (double)val[pos];
1281 means[hr][pos] += (double)val[pos];
1282 ++cell_samples[hr][pos];
1290 void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
1291 char FillFlag, float FillValue, double **means, size_t **cell_samples)
1295 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1297 for (i = 0; i < d1; i++) {
1298 #pragma _CRI concurrent
1299 for (j = 0; j < d2; j++) {
1301 if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1302 if (cell_samples[hr][pos] == 0)
1303 means[hr][pos] = val[pos];
1305 means[hr][pos] += val[pos];
1306 ++cell_samples[hr][pos];
1315 void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
1319 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1320 for (i = 0; i < d1; i++) {
1321 #pragma _CRI concurrent
1322 for (j = 0; j < d2; j++) {
1324 if (cell_samples[hr][pos] != 0) {
1325 means[hr][pos] /= (double)cell_samples[hr][pos];
1334 void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
1335 char FillFlag, float FillValue, double **means, double **stddevs,
1336 size_t **cell_samples)
1341 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1343 for (i = 0; i < d1; i++) {
1344 #pragma _CRI concurrent
1345 for (j = 0; j < d2; j++) {
1347 if (!FillFlag || (FillFlag && val[pos] != FillValue
1348 && means[hr][pos] != FillValue)) {
1349 delta = means[hr][pos] - (double)val[pos];
1350 if (cell_samples[hr][pos] == 0)
1351 stddevs[hr][pos] = delta * delta;
1353 stddevs[hr][pos] += delta * delta;
1354 ++cell_samples[hr][pos];
1362 void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
1363 char FillFlag, float FillValue, double **means, double **stddevs,
1364 size_t **cell_samples)
1369 hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1371 for (i = 0; i < d1; i++) {
1372 #pragma _CRI concurrent
1373 for (j = 0; j < d2; j++) {
1375 if (!FillFlag || (FillFlag && val[pos] != FillValue
1376 && means[hr][pos] != FillValue)) {
1377 delta = means[hr][pos] - val[pos];
1378 if (cell_samples[hr][pos] == 0)
1379 stddevs[hr][pos] = delta * delta;
1381 stddevs[hr][pos] += delta * delta;
1382 ++cell_samples[hr][pos];
1391 void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
1395 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1396 for (i = 0; i < d1; i++) {
1397 #pragma _CRI concurrent
1398 for (j = 0; j < d2; j++) {
1400 if (cell_samples[hr][pos] != 0) {
1401 stddevs[hr][pos] /= (double)cell_samples[hr][pos];
1402 stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
1411 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)
1415 size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1416 char name[NC_MAX_NAME+1];
1419 /* Allocate space for timeslice */
1420 val = alloc_var(nctype, (d1 * d2));
1422 for (i = 0; i < nifnames; i++) {
1424 printf("\tOpening %s", input_fnames[i]);
1425 if (ndims > 3) printf(" to retrieve level %d", level);
1426 if (stat_type == MEAN_TYPE)
1427 printf(" for computing means\n");
1429 printf(" for computing stddevs\n");
1431 /* Open input file */
1432 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1434 * Inquire about number of dimensions, variables, global
1435 * attributes and the ID of the unlimited dimension
1437 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1438 &input_ngatts, &input_unlimdimid));
1439 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1440 if (strcmp(name, time_name)) {
1441 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1444 /* Get seconds of day variable */
1445 wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
1446 if (!(mcsec = (int *)malloc(len * sizeof(int)))) {
1447 perror("read_and_compute:mcsec");
1450 wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));
1451 /* Get variable ID for requested variable */
1452 wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
1453 /* Retrieve time slice of desired variable */
1454 for (ts = 0; ts < len; ts++) {
1457 start[(ndims-2)] = 0;
1458 count[(ndims-2)] = d1;
1459 start[(ndims-1)] = 0;
1460 count[(ndims-1)] = d2;
1467 wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));
1468 if (stat_type == MEAN_TYPE)
1469 add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1471 add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1474 wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));
1475 if (stat_type == MEAN_TYPE)
1476 add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1478 add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1481 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1486 /* Free mcsec vector */
1488 /* Close input file */
1489 wrap_nc(nc_close(input_ncid));
1491 /* Divide sums by number of samples */
1492 if (stat_type == MEAN_TYPE)
1493 divide_means(d1, d2, means, cell_samples);
1495 divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
1497 /* Free working variable array */
1503 float *double_to_float(size_t len, double *dvar)
1508 if (!(fvar = (float *)malloc(len * sizeof(float)))) {
1509 perror("double_to_float:fvar");
1513 for (i = 0; i < len; i++)
1514 fvar[i] = (float)dvar[i];
1519 void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
1520 int level, int ndims, double **var_hours)
1522 /* Output dimensions should be one larger than input dimensions */
1524 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1527 if (nctype == NC_FLOAT) {
1528 if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
1529 perror("write_var_hours:fvar");
1534 for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1536 count[0] = 1; /* time */
1537 start[(ndims-2)] = 0;
1538 count[(ndims-2)] = d1;
1539 start[(ndims-1)] = 0;
1540 count[(ndims-1)] = d2;
1547 for (i = 0; i < (d1 * d2); i++)
1548 fvar[i] = (float)var_hours[hr][i];
1549 wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
1552 wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
1555 fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1560 if (nctype == NC_FLOAT)
1566 void compute_stats(int nifnames, char **input_fnames)
1569 size_t d1, d2, **cell_samples;
1570 double **means, **stddevs;
1571 struct var *in_vnode, *out_vnode;
1573 if (input_var_head) {
1574 in_vnode = input_var_head;
1575 /* March through non-metadata variables to compute stats */
1576 for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
1577 if (!in_vnode->metadata) {
1578 /* Make sure time is really the first dimension */
1579 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
1580 fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
1583 /* Figure out input dimensions */
1584 if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
1585 fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
1588 /* Assume 2D output; find dimensions */
1589 d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
1590 d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
1591 if (in_vnode->ndims == 3)
1594 nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
1595 /* Allocate working space for means and stddevs */
1596 alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
1597 init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
1598 printf("Computing statistics for %s\n",
1600 /* Compute means and stddevs, then write them */
1601 for (k = 0; k < nlevels; k++) {
1602 /* Read and compute means */
1603 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);
1604 /* Find corresponding output variable on the mean output file */
1605 out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
1606 /* Write out the means for this variable */
1607 write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
1608 /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
1609 reset_cell_samples(d1, d2, cell_samples);
1610 /* Read and compute stddevs using means */
1611 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);
1612 /* Find corresponding output variable on the stddev output file */
1613 out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
1614 /* Write out stddevs for this variable */
1615 write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
1618 /* Free working space for means and stddevs */
1619 free_means_stddevs(means, stddevs, cell_samples);
1621 in_vnode = in_vnode->next;
1627 void usage(char *program)
1629 fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
1630 fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
1634 int main(int argc, char **argv)
1636 int i, ic, nifnames;
1637 size_t len, pos, nsamples;
1638 char *mfname, *sfname, **ifnames, *flist;
1641 mfname = sfname = flist = NULL;
1645 printf("Number of metadata variables (nmvars) = %d\n", nmvars);
1649 /* check command-line flags and switches */
1650 while ((ic = getopt(argc, argv, "m:s:")) != -1) {
1653 if (!(mfname = strdup(optarg))) {
1659 if (!(sfname = strdup(optarg))) {
1668 fprintf(stderr, "Output file name for writing means required.\n");
1672 fprintf(stderr, "Output file name for writing standard deviations required.\n");
1676 if (optind < argc) {
1677 if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
1681 for (i = optind; i < argc; i++) {
1682 if (!(ifnames[nifnames++] = strdup(argv[i]))) {
1683 perror("ifnames[nifnames++]");
1689 fprintf(stderr, "At least one input file name is required.\n");
1695 * Build list of input files to be included in the global history
1696 * attribute on the two outputs files.
1698 for (i = len = 0; i < nifnames; i++)
1699 len += strlen(ifnames[i]);
1700 len += nifnames + 1; /* space in front of every name + null terminator */
1701 if (!(flist = (char *)malloc(len * sizeof(char)))) {
1705 for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
1706 sprintf(&flist[pos], " %s", ifnames[i]);
1708 printf("flist=%s\n", flist);
1712 /* Open every file to count up the number of time samples. */
1713 nsamples = count_samples(nifnames, ifnames);
1715 printf("Number of samples across all files = %ld\n", (long)nsamples);
1720 * Open the *last* input file and the two output files (for mean and
1721 * standard deviation). Define dimensions, variables, and attributes
1722 * for the two output files based on the *last* input files. The
1723 * last file is used because some metadata variables will contain
1724 * only values from the last time slice from the period over which
1725 * the statistics are computed.
1727 open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
1730 compute_stats(nifnames, ifnames);
1732 /* Close the two output files */
1733 wrap_nc(nc_close(mean_ncid));
1734 wrap_nc(nc_close(stddev_ncid));