forrest@0: #include <stdio.h>
forrest@0: #include <unistd.h>
forrest@0: #include <stdlib.h>
forrest@0: #include <string.h>
forrest@0: #include <math.h>
forrest@0: #include "netcdf.h"
forrest@0: 
forrest@0: /*
forrest@0:  * Loop through one month of hourly history tapes from the Community Land Model
forrest@0:  * and generate monthly statistics (means and standard deviations) of fields
forrest@0:  * by hour of day.
forrest@0:  */
forrest@0: 
forrest@0: #define	MEAN_TYPE	0
forrest@0: #define STDDEV_TYPE	1
forrest@0: 
forrest@0: #define HOURS_PER_DAY	24
forrest@0: #define MIN_PER_HOUR	60
forrest@0: #define SEC_PER_MIN	60
forrest@0: #define SEC_PER_HOUR	(MIN_PER_HOUR * SEC_PER_MIN)
forrest@0: 
forrest@0: static char *metadata_vars[] = {
forrest@0: 	"lon",
forrest@0: 	"lat",
forrest@0: 	"lonrof",
forrest@0: 	"latrof",
forrest@0: 	"time",
forrest@0: 	"hour", /* new metadata variable to be added to output files */
forrest@0: 	"levsoi",
forrest@0: 	"levlak",
forrest@0: 	"edgen",
forrest@0: 	"edgee",
forrest@0: 	"edges",
forrest@0: 	"edgew",
forrest@0: 	"longxy",
forrest@0: 	"latixy",
forrest@0: 	"area",
forrest@0: 	"areaupsc",
forrest@0: 	"topo",
forrest@0: 	"topodnsc",
forrest@0: 	"landfrac",
forrest@0: 	"landmask",
forrest@0: 	"pftmask",
forrest@0: 	"indxupsc",
forrest@0: 	"mcdate",
forrest@0: 	"mcsec",
forrest@0: 	"mdcur",
forrest@0: 	"mscur",
forrest@0: 	"nstep",
forrest@0: 	"time_bounds",
forrest@0: 	"date_written",
forrest@0: 	"time_written"
forrest@0: };
forrest@0: 
forrest@0: struct text_att {
forrest@0: 	char *name;
forrest@0: 	char *value;
forrest@0: 	struct text_att *next;
forrest@0: 	struct text_att *prev;
forrest@0: };
forrest@0: 
forrest@0: struct dim {
forrest@0: 	int dimid;
forrest@0: 	char *name;
forrest@0: 	size_t len;
forrest@0: 	struct dim *next;
forrest@0: 	struct dim *prev;
forrest@0: };
forrest@0: 
forrest@0: struct var {
forrest@0: 	int ncvarid;
forrest@0: 	char *name;
forrest@0: 	nc_type nctype;
forrest@0: 	int ndims;
forrest@0: 	int *dimids;
forrest@0: 	int natts;
forrest@0: 	char metadata;
forrest@0: 	char FillFlag;
forrest@0: 	float FillValue;
forrest@0: 	struct var *next;
forrest@0: 	struct var *prev;
forrest@0: };
forrest@0: 
forrest@0: static char *time_name = "time";
forrest@1: static char *time_bounds_name = "time_bounds";
forrest@0: static char *mcsec_name = "mcsec";
forrest@0: static char *history_name = "history";
forrest@0: static char *nsamples_name = "nsamples";
forrest@0: static char *hour_name = "hour", *hour_long_name = "hour of day",
forrest@0: 	*hour_units = "hour";
forrest@0: static char *cell_method_name = "cell_method";
forrest@0: /* input stuff */
forrest@0: static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
forrest@0: static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
forrest@0: static struct text_att *input_gatt_head = NULL;
forrest@0: static struct dim *input_dim_head = NULL, **input_dim_idx;
forrest@0: static struct var *input_var_head = NULL;
forrest@0: /* translation stuff */
forrest@0: static int *idid2mdid, *idid2sdid; /* dimension IDs */
forrest@0: /* output stuff */
forrest@0: static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
forrest@0: static int mean_hour_dimid; /* special notes */
forrest@0: static struct dim *mean_dim_head = NULL;
forrest@0: static struct var *mean_var_head = NULL;
forrest@0: static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
forrest@0: static int stddev_hour_dimid; /* special notes */
forrest@0: static struct dim *stddev_dim_head = NULL;
forrest@0: static struct var *stddev_var_head = NULL;
forrest@0: 
forrest@0: char is_metadata(char *name)
forrest@0: {
forrest@0: 	int i;
forrest@0: 
forrest@0: 	for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
forrest@0: 
forrest@0: 	if (i < nmvars)
forrest@0: 		return 1;
forrest@0: 	else
forrest@0: 		return 0;
forrest@0: }
forrest@0: 
forrest@0: #if 0
forrest@0: struct dim *dimlist_find_by_name(struct dim *head, char *name)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	struct dim *p = NULL;
forrest@0: 
forrest@0: 	if (head) {
forrest@0: 		node = head;
forrest@0: 		for (i = 0 ; i == 0 || node != head; i++) {
forrest@0: 			if (!strcmp(node->name, name))
forrest@0: 				p = node;
forrest@0: 			node = node->next;
forrest@0: 		}
forrest@0: 		return p;
forrest@0: 	}
forrest@0: 	else
forrest@0: 		return NULL;
forrest@0: }
forrest@0: #endif
forrest@0: 
forrest@0: struct var *varlist_find_by_name(struct var *head, char *name)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	struct var *node;
forrest@0: 
forrest@0: 	if (head) {
forrest@0: 		node = head;
forrest@0: 		for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
forrest@0: 			node = node->next;
forrest@0: 		if (i && node == head)
forrest@0: 			return NULL;
forrest@0: 		else
forrest@0: 			return node;
forrest@0: 	}
forrest@0: 	else
forrest@0: 		return NULL;
forrest@0: }
forrest@0: 
forrest@0: void gattlist_add_node(struct text_att **headp, char *name, char *value)
forrest@0: {
forrest@0: 	struct text_att *head, *node;
forrest@0: 
forrest@0: 	head = *headp;
forrest@0: 
forrest@0: 	if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
forrest@0: 		perror("gattlist_add_node");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	if (!(node->name = strdup(name))) {
forrest@0: 		perror("gattlist_add_node");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	if (!(node->value = strdup(value))) {
forrest@0: 		perror("gattlist_add_node");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	if (head == NULL) {
forrest@0: 		/* first one */
forrest@0: 		*headp = node;
forrest@0: 		node->next = node;
forrest@0: 		node->prev = node;
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		node->next = head;
forrest@0: 		node->prev = head->prev;
forrest@0: 		/* set this after setting node->prev from here */
forrest@0: 		head->prev = node;
forrest@0: 		/* set this after having set node->prev */
forrest@0: 		node->prev->next = node;
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
forrest@0: {
forrest@0: 	struct dim *head, *node;
forrest@0: 
forrest@0: 	head = *headp;
forrest@0: 
forrest@0: 	if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
forrest@0: 		perror("dimlist_add_node");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	node->dimid = dimid;
forrest@0: 	if (!(node->name = strdup(name))) {
forrest@0: 		perror("dimlist_add_node");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	node->len = len;
forrest@0: 	
forrest@0: 	if (head == NULL) {
forrest@0: 		/* first one */
forrest@0: 		*headp = node;
forrest@0: 		node->next = node;
forrest@0: 		node->prev = node;
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		node->next = head;
forrest@0: 		node->prev = head->prev;
forrest@0: 		/* set this after setting node->prev from here */
forrest@0: 		head->prev = node;
forrest@0: 		/* set this after having set node->prev */
forrest@0: 		node->prev->next = node;
forrest@0: 	}
forrest@0: 
forrest@0: 	return node;
forrest@0: }
forrest@0: 
forrest@0: void varlist_add_node(struct var **headp, int ncvarid, char *name,
forrest@0: 	nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
forrest@0: 	float FillValue)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	struct var *head, *node;
forrest@0: 
forrest@0: 	head = *headp;
forrest@0: 
forrest@0: 	if (!(node = (struct var *)malloc(sizeof(struct var)))) {
forrest@0: 		perror("varlist_add_node");
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	node->ncvarid = ncvarid;
forrest@0: 	if (!(node->name = strdup(name))) {
forrest@0: 		perror("varlist_add_node");
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 	node->nctype = nctype;
forrest@0: 	node->ndims = ndims;
forrest@0: 	if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
forrest@0: 		perror("varlist_add_node");
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 	for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
forrest@0: 	node->natts = natts;
forrest@0: 	node->metadata = is_metadata(name);
forrest@0: 	node->FillFlag = FillFlag;
forrest@0: 	node->FillValue = FillValue;
forrest@0: 	
forrest@0: 	if (head == NULL) {
forrest@0: 		/* first one */
forrest@0: 		*headp = node;
forrest@0: 		node->next = node;
forrest@0: 		node->prev = node;
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		node->next = head;
forrest@0: 		node->prev = head->prev;
forrest@0: 		/* set this after setting node->prev from here */
forrest@0: 		head->prev = node;
forrest@0: 		/* set this after having set node->prev */
forrest@0: 		node->prev->next = node;
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void varlist_print(struct var *headp)
forrest@0: {
forrest@0: 	int i, j;
forrest@0: 	struct var *node;
forrest@0: 
forrest@0: 	printf("Beginning of Variable List\n");
forrest@0: 
forrest@0: 	if (headp) {
forrest@0: 		node = headp;
forrest@0: 		for (j = 0; j == 0 || node != headp; j++) {
forrest@0: 			printf("Variable %d (ptr=%ld):\n", j, (long)node);
forrest@0: 			printf("\tncvarid = %d\n", node->ncvarid);
forrest@0: 			printf("\tname = %s\n",  node->name);
forrest@0: 			printf("\tnctype = %d\n", node->nctype);
forrest@0: 			printf("\tndims = %d\n", node->ndims);
forrest@0: 			printf("\tdimids =");
forrest@0: 			for (i = 0; i < node->ndims; i++)
forrest@0: 				printf(" %d", node->dimids[i]);
forrest@0: 			printf("\n");
forrest@0: 			printf("\tnatts = %d\n", node->natts);
forrest@0: 			printf("\tmetadata = %d\n", (int)node->metadata);
forrest@0: 			printf("\tnext = %ld\n", (long)node->next);
forrest@0: 			printf("\tprev = %ld\n", (long)node->prev);
forrest@0: 			node = node->next;
forrest@0: 		}
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		printf("\tList undefined\n");
forrest@0: 	}
forrest@0: 
forrest@0: 	printf("End of Variable List\n");
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void wrap_nc(int status)
forrest@0: {
forrest@0: 	if (status != NC_NOERR) {
forrest@0: 		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
forrest@0: 		exit(1);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	char name[NC_MAX_NAME+1];
forrest@0: 	size_t len;
forrest@0: 	struct dim **in_dim_idx;
forrest@0: 
forrest@0: 	if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
forrest@0: 		perror("*in_dim_idxp");
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 	in_dim_idx = *in_dim_idxp;
forrest@0: 
forrest@0: 	for (i = 0; i < ndims; i++) {
forrest@0: 		wrap_nc(nc_inq_dim(ncid, i, name, &len));
forrest@0: 		in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	char name[NC_MAX_NAME+1], *value;
forrest@0: 	nc_type xtype;
forrest@0: 	size_t len;
forrest@0: 
forrest@0: 	for (i = 0; i < ngatts; i++) {
forrest@0: 		wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
forrest@0: 		wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
forrest@0: 		if (xtype != NC_CHAR) {
forrest@0: 			fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 		if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
forrest@0: 			perror("get_gatts: value");
forrest@0: 			exit(3);
forrest@0: 		}
forrest@0: 		wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
forrest@0: 		value[(len+1-1)] = '\0';
forrest@0: 		gattlist_add_node(gatt_headp, name, value);
forrest@0: 		free(value); /* because gattlist_add_node() duplicates it */
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void get_vars(int ncid, int nvars, struct var **var_headp)
forrest@0: {
forrest@0: 	int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
forrest@0: 	size_t len;
forrest@0: 	float FillValue;
forrest@0: 	char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
forrest@0: 	nc_type xtype, atype;
forrest@0: 
forrest@0: 	for (i = 0; i < nvars; i++) {
forrest@0: 		wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
forrest@0: 			&natts));
forrest@0: 		/* Check for _FillValue */
forrest@0: 		FillFlag = 0;
forrest@0: 		FillValue = 0.;
forrest@0: 		if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
forrest@0: 			== NC_NOERR) {
forrest@0: 			if (atype == NC_FLOAT && len) {
forrest@0: 				wrap_nc(nc_get_att_float(ncid, i,
forrest@0: 					fill_att_name, &FillValue));
forrest@0: 				FillFlag = 1;
forrest@0: 			}
forrest@0: 		}
forrest@0: 		varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
forrest@0: 			natts, FillFlag, FillValue);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
forrest@0: 	size_t nsamples, int **in2outp, int out_ncid,
forrest@0: 	struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Define dimensions on new files and build translation tables between
forrest@0: 	 * dimension IDs.
forrest@0: 	 */
forrest@0: 	int j, dimid, ndims, *in2out;
forrest@0: 	size_t len;
forrest@0: 	struct dim *dnode;
forrest@0: 
forrest@0: 	ndims = 0;
forrest@0: 
forrest@0: 	if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
forrest@0: 		perror("put_dimensions: in2outp");
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 	in2out = *in2outp;
forrest@0: 
forrest@0: 	if (in_dim_head) {
forrest@0: 		dnode = in_dim_head;
forrest@0: 		for (j = 0; j == 0 || dnode != in_dim_head; j++) {
forrest@0: 			if (dnode->dimid == in_unlimdimid)
forrest@0: 				len = NC_UNLIMITED;
forrest@0: 			else
forrest@0: 				len = dnode->len;
forrest@0: 			wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
forrest@0: 			dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
forrest@0: 			++ndims;
forrest@0: 			in2out[dnode->dimid] = dimid;
forrest@0: 			if (dnode->dimid == in_unlimdimid)
forrest@0: 				*out_unlimdimidp = dimid;
forrest@0: 			/*
forrest@0: 			 * Just after the "time" dimension, add the new
forrest@0: 			 * "nsamples" and "hour" dimensions.
forrest@0: 			 */
forrest@0: 			if (!strcmp(dnode->name, time_name)) {
forrest@0: 				wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
forrest@0: 				dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
forrest@0: 				++ndims;
forrest@0: 
forrest@0: 				wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
forrest@0: 				dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
forrest@0: 				++ndims;
forrest@0: 				/* Track hour dimid for out files */
forrest@0: 				*out_hour_dimidp = dimid;
forrest@0: 			}
forrest@0: 
forrest@0: 			dnode = dnode->next;
forrest@0: 		}
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		fprintf(stderr, "WARNING: No dimensions defined!\n");
forrest@0: 	}
forrest@0: 
forrest@0: 	return ndims;
forrest@0: }
forrest@0: 
forrest@0: int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Write out global attributes matching those of the input file.
forrest@0: 	 * Change history attribute to the string passed in as an argument.
forrest@0: 	 */
forrest@0: 	int j, hflag = 0, ngatts;
forrest@0: 	char *value;
forrest@0: 	struct text_att *anode;
forrest@0: 
forrest@0: 	ngatts = 0;
forrest@0: 
forrest@0: 	if (in_gatt_head) {
forrest@0: 		anode = in_gatt_head;
forrest@0: 		for (j = 0; j == 0 || anode != in_gatt_head; j++) {
forrest@0: 			if (!strcmp(anode->name, history_name)) {
forrest@0: 				value = out_history;
forrest@0: 				++hflag;
forrest@0: 			}
forrest@0: 			else
forrest@0: 				value = anode->value;
forrest@0: 			wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
forrest@0: 			++ngatts;
forrest@0: 			anode = anode->next;
forrest@0: 		}
forrest@0: 		/* If no history attribute on input, add one to the output */
forrest@0: 		if (!hflag) {
forrest@0: 			wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
forrest@0: 			++ngatts;
forrest@0: 		}
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		fprintf(stderr, "WARNING: No global attributes defined!\n");
forrest@0: 	}
forrest@0: 
forrest@0: 	return ngatts;
forrest@0: }
forrest@0: 
forrest@0: int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Translate between input dimension IDs and those for the output file.
forrest@0: 	 * For normal time-based variables, add a new dimension for hour of
forrest@0: 	 * day.  For metadata variables, do not add new dimensions, even if
forrest@0: 	 * it is time-based.  Return (possibly new) number of dimensions.
forrest@0: 	 */
forrest@0: 	int i, ndims;
forrest@0: 
forrest@0: 	for (i = ndims = 0; i < in_ndims; i++, ndims++) {
forrest@0: 		out_dimids[ndims] = in2out[in_dimids[i]];
forrest@0: 		if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
forrest@0: 			&& !metadata) {
forrest@0: 			ndims++;
forrest@0: 			out_dimids[ndims] = hour_dimid;
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return ndims;
forrest@0: }
forrest@0: 
forrest@0: int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
forrest@0: 	int in_varid, int out_ncid, int out_varid)
forrest@0: {
forrest@0: 	/* 
forrest@0: 	 * Copy the attributes of the input variable to those of the output
forrest@0: 	 * variable. If the variable is not a metadata variable, ensure that
forrest@0: 	 * the cell_method attribute is set either to "time: mean" or
forrest@0: 	 * "time: standard_deviation", depending on the output file type.
forrest@0: 	 */
forrest@0: 
forrest@0: 	int i, natts;
forrest@0: 	char name[NC_MAX_NAME+1], cmflag = 0;
forrest@0: 	char *cell_methods[] = { "time: mean", "time: standard_deviation" };
forrest@0: 
forrest@0: 	for (i = natts = 0; i < in_natts; i++, natts++) {
forrest@0: 		wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
forrest@0: 		if (!strcmp(name, cell_method_name) && !metadata) {
forrest@0: 			wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
forrest@0: 			cmflag = 1;
forrest@0: 		}
forrest@0: 		else
forrest@0: 			wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
forrest@0: 	}
forrest@0: 	/*
forrest@0: 	 * If no cell_method attribute was specified for a non-metadata
forrest@0: 	 * variable on the input file, add the appropriate cell_method anyway
forrest@0: 	 * on the output file.
forrest@0: 	 */
forrest@0: 	if (!cmflag && !metadata) {
forrest@0: 		wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
forrest@0: 		natts++;
forrest@0: 	}
forrest@0: 
forrest@0: 	return natts;
forrest@0: }
forrest@0: 
forrest@0: int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
forrest@0: 	int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp,
forrest@0: 	int stat_type)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Define variables on output file and copy attributes from input file.
forrest@0: 	 * Return number of variables defined.
forrest@0: 	 */
forrest@0: 	int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
forrest@0: 	struct var *vnode;
forrest@0: 
forrest@0: 	nvars = 0;
forrest@0: 
forrest@0: 	if (in_var_head) {
forrest@0: 		vnode = in_var_head;
forrest@0: 		/*
forrest@0: 		 * March through input variables creating (mostly) the same
forrest@0: 		 * variables on the output file.
forrest@0: 		 */
forrest@0: 		for (j = 0; j == 0 || vnode != in_var_head; j++) {
forrest@0: 			ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
forrest@0: 			wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
forrest@0: 			natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
forrest@0: 			varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
forrest@0: 			++nvars;
forrest@0: 			/*
forrest@0: 			 * Just after the "time" variable, add the new "hour"
forrest@0: 			 * variable for hour of day.
forrest@0: 			 */
forrest@0: 			if (!strcmp(vnode->name, time_name)) {
forrest@0: 				ndims = 1;
forrest@0: 				dimids[0] = hour_dimid;
forrest@0: 				wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
forrest@0: 				wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
forrest@0: 				wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
forrest@0: 				varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
forrest@0: 				++nvars;
forrest@0: 			}
forrest@0: 
forrest@0: 			vnode = vnode->next;
forrest@0: 		}
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		fprintf(stderr, "WARNING: No variables defined!\n");
forrest@0: 	}
forrest@0: 
forrest@0: 	return nvars;
forrest@0: }
forrest@0: 
forrest@0: void *alloc_var(nc_type nctype, size_t len)
forrest@0: {
forrest@0: 	void *val;
forrest@0: 
forrest@0: 	switch (nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			if (!(val = (float *)malloc((len) * sizeof(float)))) {
forrest@0: 				perror("alloc_var: val");
forrest@0: 				exit(3);
forrest@0: 			}
forrest@0: 			break;
forrest@0: 		case NC_INT:
forrest@0: 			if (!(val = (int *)malloc((len) * sizeof(int)))) {
forrest@0: 				perror("alloc_var: val");
forrest@0: 				exit(3);
forrest@0: 			}
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			if (!(val = (double *)malloc((len) * sizeof(double)))) {
forrest@0: 				perror("alloc_var: val");
forrest@0: 				exit(3);
forrest@0: 			}
forrest@0: 			break;
forrest@0: 		case NC_CHAR:
forrest@0: 			if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
forrest@0: 				perror("alloc_var: val");
forrest@0: 				exit(3);
forrest@0: 			}
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	return val;
forrest@0: }
forrest@0: 
forrest@0: void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	size_t len = (size_t)1;
forrest@0: 	void *val;
forrest@0: 
forrest@0: 	/* Figure out the total size */
forrest@0: 	for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
forrest@0: 
forrest@0: 	val = alloc_var(nctype, len);
forrest@0: 
forrest@0: 	switch (nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			wrap_nc(nc_get_var_float(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_INT:
forrest@0: 			wrap_nc(nc_get_var_int(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			wrap_nc(nc_get_var_double(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_CHAR:
forrest@0: 			wrap_nc(nc_get_var_text(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	return val;
forrest@0: }
forrest@0: 
forrest@0: void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 	void *val;
forrest@0: 
forrest@0: 	/* Make sure time is really the first dimension */
forrest@0: 	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
forrest@0: 		fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
forrest@0: 		exit(3);
forrest@0: 	} 
forrest@0: 	/*
forrest@0: 	 * Figure out the total size for the slice (start at second dimension)
forrest@0: 	 * and build start/count vectors.
forrest@0: 	 */
forrest@0: 	start[0] = tslice;
forrest@0: 	count[0] = 1;
forrest@0: 	for (i = 1; i < ndims; i++) {
forrest@0: 		start[i] = 0;
forrest@0: 		count[i] = (dim_idx[dimids[i]])->len;
forrest@0: 		len *= (dim_idx[dimids[i]])->len;
forrest@0: 	}
forrest@0: 
forrest@0: 	val = alloc_var(nctype, len);
forrest@0: 	
forrest@0: 	switch (nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_INT:
forrest@0: 			wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_CHAR:
forrest@0: 			wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	return val;
forrest@0: }
forrest@0: 
forrest@0: void write_var(int ncid, int varid, nc_type nctype, void *val)
forrest@0: {
forrest@0: 	switch (nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			wrap_nc(nc_put_var_float(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_INT:
forrest@0: 			wrap_nc(nc_put_var_int(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			wrap_nc(nc_put_var_double(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		case NC_CHAR:
forrest@0: 			wrap_nc(nc_put_var_text(ncid, varid, val));
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 
forrest@0: 	/* Make sure time is really the first dimension */
forrest@0: 	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
forrest@0: 		fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
forrest@0: 		exit(3);
forrest@0: 	}
forrest@0: 	
forrest@0: 	/* Build start/count vectors */
forrest@0: 	start[0] = tslice;
forrest@0: 	count[0] = 1;
forrest@0: 	for (i = 1; i < ndims; i++) {
forrest@0: 		start[i] = 0;
forrest@0: 		count[i] = (dim_idx[dimids[i]])->len;
forrest@0: 	}
forrest@0: 
forrest@0: 	switch (nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_INT:
forrest@0: 			wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		case NC_CHAR:
forrest@0: 			wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void copy_metadata(int in_ncid, struct var *in_var_head,
forrest@1: 	struct dim **in_dim_idx, int out_ncid, struct var *out_var_head,
forrest@1: 	double *tbounds)
forrest@0: {
forrest@0: 	int i, j;
forrest@0: 	float hr[HOURS_PER_DAY];
forrest@0: 	struct var *in_vnode, *out_vnode;
forrest@0: 	void *val;
forrest@0: 
forrest@0: 	for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
forrest@0: 
forrest@0: 	if (in_var_head) {
forrest@0: 		in_vnode = in_var_head;
forrest@0: 		/*
forrest@0: 		 * March through input variables to stuff metadata values into
forrest@0: 		 * the new output files. NOTE: Time-based metadata variables
forrest@0: 		 * should contain only the last (time-slice) value (from all
forrest@0: 		 * input files).
forrest@0: 		 */
forrest@0: 		for (j = 0; j == 0 || in_vnode != in_var_head; j++) {
forrest@0: 			if (in_vnode->metadata) {
forrest@0: 				out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
forrest@0: 				if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
forrest@0: 					/* time is not the first dimension */
forrest@0: #ifdef DEBUG
forrest@0: 					printf("Copying metadata variable %s\n",
forrest@0: 						in_vnode->name);
forrest@0: #endif /* DEBUG */
forrest@0: 					val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
forrest@0: 					write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
forrest@0: 				}
forrest@0: 				else {
forrest@0: 					/* time is the first dimension */
forrest@0: #ifdef DEBUG
forrest@0: 					printf("Copying last value of \
forrest@0: time-based metadata variable %s\n", in_vnode->name);
forrest@0: #endif /* DEBUG */
forrest@0: 					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));
forrest@1: 					if (!strcmp(in_vnode->name, time_bounds_name))
forrest@1: 						write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, tbounds, 0);
forrest@1: 					else
forrest@1: 						write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
forrest@0: 				}
forrest@0: 				free(val);
forrest@0: 				/*
forrest@0: 				 * Just after the "time" variable, write out
forrest@0: 				 * the "hour" variable values.
forrest@0: 				 */
forrest@0: 				if (!strcmp(in_vnode->name, time_name)) {
forrest@0: 					out_vnode = varlist_find_by_name(out_var_head, hour_name);
forrest@0: 					write_var(out_ncid, out_vnode->ncvarid,
forrest@0: 						out_vnode->nctype, hr);
forrest@0: 				}
forrest@0: 			}
forrest@0: 			else {
forrest@0: #ifdef DEBUG
forrest@0: 				printf("Skipping variable %s\n",
forrest@0: 					in_vnode->name);
forrest@0: #endif /* DEBUG */
forrest@0: 			}
forrest@0: 			in_vnode = in_vnode->next;
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@1: void get_time_bounds(char *first_fname, char *last_fname, double *tbounds)
forrest@1: {
forrest@1: 	int time_dimid, time_bounds_varid;
forrest@1: 	size_t time_len, time_bounds_index[2];
forrest@1: 	double bnd1, bnd2;
forrest@1: 
forrest@1: 	/* Open first input file */
forrest@1: 	wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
forrest@1: 	/* Get dimension ID of the time dimension */
forrest@1: 	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1: 	/* Get length of time dimension */
forrest@1: 	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1: 	/* Get variable ID of time_bounds variable */
forrest@1: 	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
forrest@1: 	/* Set index for reading the first value of time_bounds from the
forrest@1: 	 * first timeslice of the first file */
forrest@1: 	time_bounds_index[0] = time_bounds_index[1] = 0;
forrest@1: 	/* Read the value */
forrest@1: 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1: 		time_bounds_index, &bnd1));
forrest@1: 	/* If the first and last file are not the same, close the first one and
forrest@1: 	 * open the second one */
forrest@1: 	if (strcmp(first_fname, last_fname)) {
forrest@1: 		/* Close the first input file */
forrest@1: 		wrap_nc(nc_close(input_ncid));
forrest@1: 		/* Open the last input file */
forrest@1: 		wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@1: 		/* Get dimension ID of the time dimension */
forrest@1: 		wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1: 		/* Get length of time dimension */
forrest@1: 		wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1: 		/* Get variable ID of time_bounds variable */
forrest@1: 		wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
forrest@1: 			&time_bounds_varid));
forrest@1: 	}
forrest@1: 	/* Set index for reading the second value of time_bounds from the
forrest@1: 	 * last timeslice of the last file */
forrest@1: 	time_bounds_index[0] = time_len - 1; time_bounds_index[1] = 1;
forrest@1: 	/* Read the value */
forrest@1: 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1: 		time_bounds_index, &bnd2));
forrest@1: 
forrest@1: 	/* Close the last input file */
forrest@1: 	wrap_nc(nc_close(input_ncid));
forrest@1: 
forrest@1: 	tbounds[0] = bnd1;
forrest@1: 	tbounds[1] = bnd2;
forrest@1: 
forrest@1: 	return;
forrest@1: }
forrest@1: 
forrest@1: void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
forrest@0: {
forrest@0: 	char *mean_history_gatt, *stddev_history_gatt,
forrest@0: 		*mean_prefix="Statistical means from history files:",
forrest@0: 		*stddev_prefix="Statistical standard deviations from history files:";
forrest@1: 	double tbounds[2];
forrest@0: 
forrest@0: 	/*
forrest@0: 	 * Construct strings for history global attributes for the two output
forrest@0: 	 * files.
forrest@0: 	 */
forrest@0: 	if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
forrest@0: 		perror("open_inout:mean_history_gatt");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
forrest@0: 		perror("open_inout:stddev_history_gatt");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
forrest@0: 	sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
forrest@0: 
forrest@1: 	/* The two time_bounds values must be handled explicitly by obtaining
forrest@1: 	 * the first one from the first file and the last one from the last
forrest@1: 	 * file.  These are then passed to copy_metadata() below. */
forrest@1: 	get_time_bounds(first_fname, last_fname, tbounds);
forrest@1: #ifdef DEBUG
forrest@1: 	fprintf(stderr, "Got back tbounds=%lf,%lf\n", tbounds[0], tbounds[1]);
forrest@1: #endif /* DEBUG */
forrest@1: 
forrest@1: 	/* Open last input file */
forrest@1: 	wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@0: 	/* Inquire about number of dimensions, variables, global attributes
forrest@0: 	 * and the ID of the unlimited dimension
forrest@0: 	 */
forrest@0: 	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
forrest@0: 		&input_unlimdimid));
forrest@0: 	/* Grab dimension IDs and lengths from input file */
forrest@0: 	get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
forrest@0: 	/* Grab desired global attributes from input file */
forrest@0: 	get_gatts(input_ncid, input_ngatts, &input_gatt_head);
forrest@0: 	/* Grab variable IDs and attributes from input file */
forrest@0: 	get_vars(input_ncid, input_nvars, &input_var_head);
forrest@0: #ifdef DEBUG
forrest@0: 	varlist_print(input_var_head);
forrest@0: #endif /* DEBUG */
forrest@0: 	/* Create netCDF files for monthly means and standard deviations */
forrest@0: 	/* Create new netCDF files */
forrest@0: 	wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
forrest@0: 	wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
forrest@0: 	/* Define dimensions */
forrest@0: 	mean_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0: 		input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
forrest@0: 		&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
forrest@0: 	stddev_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0: 		input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
forrest@0: 		&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
forrest@0: 	/* Define variables and their attributes */
forrest@0: 	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0: 		idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
forrest@0: 		MEAN_TYPE);
forrest@0: 	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0: 		idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
forrest@0: 		STDDEV_TYPE);
forrest@0: 	/* Store global attributes */
forrest@0: 	mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
forrest@0: 	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
forrest@0: 		stddev_history_gatt);
forrest@0: 	/* End define mode */
forrest@0: 	wrap_nc(nc_enddef(mean_ncid));
forrest@0: 	wrap_nc(nc_enddef(stddev_ncid));
forrest@0: 	/* Write out metdata variables that do not depend on "time" */
forrest@0: 	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
forrest@1: 		mean_var_head, tbounds);
forrest@0: 	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
forrest@1: 		stddev_var_head, tbounds);
forrest@0: 
forrest@0: 	wrap_nc(nc_close(input_ncid));
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: size_t count_samples(int nifnames, char **input_fnames)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	char name[NC_MAX_NAME+1];
forrest@0: 	size_t len, cnt = 0;
forrest@0: 
forrest@0: 	for (i = 0; i < nifnames; i++) {
forrest@0: 		/* Open input file */
forrest@0: 		wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0: 		/*
forrest@0: 		 * Inquire about number of dimensions, variables, global
forrest@0: 		 * attributes and the ID of the unlimited dimension
forrest@0: 		 */
forrest@0: 		wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0: 			&input_ngatts, &input_unlimdimid));
forrest@0: 		wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0: 		if (strcmp(name, time_name)) {
forrest@0: 			fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0: 			exit(4);
forrest@0: 		}
forrest@0: #ifdef DEBUG
forrest@0: 		printf("%ld samples in %s\n", (long)len, input_fnames[i]);
forrest@0: #endif /* DEBUG */
forrest@0: 		wrap_nc(nc_close(input_ncid));
forrest@0: 		cnt += len;
forrest@0: 	}
forrest@0: 	return cnt;
forrest@0: }
forrest@0: 
forrest@0: void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Allocate space for arrays of means and standard deviations by
forrest@0: 	 * hour of day.
forrest@0: 	 */
forrest@0: 	int i;
forrest@0: 	size_t **cell_samples;
forrest@0: 	double **means, **stddevs;
forrest@0: 
forrest@0: 	if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
forrest@0: 		perror("alloc_means_stddevs:*meansp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
forrest@0: 		perror("alloc_means_stddevs:*stddevsp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
forrest@0: 		perror("alloc_means_stddevs:*cell_samplesp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	means = *meansp;
forrest@0: 	stddevs = *stddevsp;
forrest@0: 	cell_samples = *cell_samplesp;
forrest@0: 
forrest@0: 	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0: 		if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
forrest@0: 			perror("alloc_means_stddevs:means[i]");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 		if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
forrest@0: 			perror("alloc_means_stddevs:stddevs[i]");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 		if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
forrest@0: 			perror("alloc_means_stddevs:cell_samples[i]");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
forrest@0: {
forrest@0: 	int hr, i, j, pos;
forrest@0: 
forrest@0: 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0: 		for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 			for (j = 0; j < d2; j++) {
forrest@0: 				pos = i*d2+j;
forrest@0: 				means[hr][pos] = FillValue;
forrest@0: 				stddevs[hr][pos] = FillValue;
forrest@0: 				cell_samples[hr][pos] = 0;
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 
forrest@0: 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0: 		for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 			for (j = 0; j < d2; j++) {
forrest@0: 				pos = i*d2+j;
forrest@0: 				cell_samples[hr][pos] = 0;
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
forrest@0: {
forrest@0: 	/*
forrest@0: 	 * Free space from arrays of means and standard deviations by
forrest@0: 	 * hour of day.
forrest@0: 	 */
forrest@0: 	int i;
forrest@0: 
forrest@0: 	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0: 		free(means[i]);
forrest@0: 		free(stddevs[i]);
forrest@0: 		free(cell_samples[i]);
forrest@0: 	}
forrest@0: 
forrest@0: 	free(means);
forrest@0: 	free(stddevs);
forrest@0: 	free(cell_samples);
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
forrest@0: 	char FillFlag, float FillValue, double **means, size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 
forrest@0: 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0: 
forrest@0: 	for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 		for (j = 0; j < d2; j++) {
forrest@0: 			pos = i*d2+j;
forrest@0: 			if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
forrest@0: 				if (cell_samples[hr][pos] == 0)
forrest@0: 					means[hr][pos] = (double)val[pos];
forrest@0: 				else
forrest@0: 					means[hr][pos] += (double)val[pos];
forrest@0: 				++cell_samples[hr][pos];
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
forrest@0: 	char FillFlag, float FillValue, double **means, size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 
forrest@0: 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0: 
forrest@0: 	for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 		for (j = 0; j < d2; j++) {
forrest@0: 			pos = i*d2+j;
forrest@0: 			if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
forrest@0: 				if (cell_samples[hr][pos] == 0)
forrest@0: 					means[hr][pos] = val[pos];
forrest@0: 				else
forrest@0: 					means[hr][pos] += val[pos];
forrest@0: 				++cell_samples[hr][pos];
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: 
forrest@0: void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 
forrest@0: 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0: 		for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 			for (j = 0; j < d2; j++) {
forrest@0: 				pos = i*d2+j;
forrest@0: 				if (cell_samples[hr][pos] != 0) {
forrest@0: 					means[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0: 				}
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
forrest@0: 	char FillFlag, float FillValue, double **means, double **stddevs,
forrest@0: 	size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 	double delta;
forrest@0: 
forrest@0: 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0: 
forrest@0: 	for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 		for (j = 0; j < d2; j++) {
forrest@0: 			pos = i*d2+j;
forrest@0: 			if (!FillFlag || (FillFlag && val[pos] != FillValue
forrest@0: 				&& means[hr][pos] != FillValue)) {
forrest@0: 				delta = means[hr][pos] - (double)val[pos];
forrest@0: 				if (cell_samples[hr][pos] == 0)
forrest@0: 					stddevs[hr][pos] = delta * delta;
forrest@0: 				else
forrest@0: 					stddevs[hr][pos] += delta * delta;
forrest@0: 				++cell_samples[hr][pos];
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
forrest@0: 	char FillFlag, float FillValue, double **means, double **stddevs,
forrest@0: 	size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 	double delta;
forrest@0: 
forrest@0: 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0: 
forrest@0: 	for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 		for (j = 0; j < d2; j++) {
forrest@0: 			pos = i*d2+j;
forrest@0: 			if (!FillFlag || (FillFlag && val[pos] != FillValue
forrest@0: 				&& means[hr][pos] != FillValue)) {
forrest@0: 				delta = means[hr][pos] - val[pos];
forrest@0: 				if (cell_samples[hr][pos] == 0)
forrest@0: 					stddevs[hr][pos] = delta * delta;
forrest@0: 				else
forrest@0: 					stddevs[hr][pos] += delta * delta;
forrest@0: 				++cell_samples[hr][pos];
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: 
forrest@0: void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
forrest@0: {
forrest@0: 	int i, j, hr, pos;
forrest@0: 
forrest@0: 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0: 		for (i = 0; i < d1; i++) {
forrest@0: #pragma _CRI concurrent
forrest@0: 			for (j = 0; j < d2; j++) {
forrest@0: 				pos = i*d2+j;
forrest@0: 				if (cell_samples[hr][pos] != 0) {
forrest@0: 					stddevs[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0: 					stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
forrest@0: 				}
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	void **vals;
forrest@0: 
forrest@0: 	/* Allocate space for all timeslices */
forrest@0: 	if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) {
forrest@0: 		perror("alloc_all_samples:*valsp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	vals = *valsp;
forrest@0: 	for (i = 0; i < nsamples; i++) {
forrest@0: 		vals[i] = alloc_var(nctype, (d1 * d2));
forrest@0: 	}
forrest@0: 	if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) {
forrest@0: 		perror("alloc_all_samples:*mcsecp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec)
forrest@0: {
forrest@0: 	int i, j;
forrest@0:         float **fvals = NULL;
forrest@0: 	double **dvals = NULL;
forrest@0: 
forrest@0: 	switch(nctype) {
forrest@0: 		case NC_FLOAT:
forrest@0: 			fvals = (float **)vals;
forrest@0: 			break;
forrest@0: 		case NC_DOUBLE:
forrest@0: 			dvals = (double **)vals;
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 
forrest@0: 	for (i = 0; i < nsamples; i++) {
forrest@0: 		for (j = 0; j < (d1 * d2); j++) {
forrest@0: 			switch(nctype) {
forrest@0: 				case NC_FLOAT:
forrest@0: 					fvals[i][j] = (float)0.0;
forrest@0: 					break;
forrest@0: 				case NC_DOUBLE:
forrest@0: 					dvals[i][j] = (double)0.0;
forrest@0: 					break;
forrest@0: 				default:
forrest@0: 					fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 					exit(3);
forrest@0: 			}
forrest@0: 		}
forrest@0: 		mcsec[i] = 0;
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void free_all_samples(size_t nsamples, void **vals, int *mcsec)
forrest@0: {
forrest@0: 	int i;
forrest@0: 
forrest@0: 	for (i = 0; i < nsamples; i++)
forrest@0: 		if (vals[i]) free(vals[i]);
forrest@0: 
forrest@0: 	free(vals);
forrest@0: 	free(mcsec);
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: 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)
forrest@0: {
forrest@0: 	int i, ts;
forrest@0: 	int varid;
forrest@0: 	size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 	char name[NC_MAX_NAME+1];
forrest@0: 
forrest@0: 	for (i = 0; i < nifnames; i++) {
forrest@0: #ifdef DEBUG
forrest@0: 		printf("\tOpening %s", input_fnames[i]);
forrest@0: 		if (ndims > 3) printf(" to retrieve level %d\n", level);
forrest@0: #endif /* DEBUG */
forrest@0: 		/* Open input file */
forrest@0: 		wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0: 		/*
forrest@0: 		 * Inquire about number of dimensions, variables, global
forrest@0: 		 * attributes and the ID of the unlimited dimension
forrest@0: 		 */
forrest@0: 		wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0: 			&input_ngatts, &input_unlimdimid));
forrest@0: 		wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0: 		if (strcmp(name, time_name)) {
forrest@0: 			fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0: 			exit(4);
forrest@0: 		}
forrest@0: 		/* Make sure we don't run off the end of allocated space */
forrest@0: 		if ((nslice+len) > nsamples) {
forrest@0: 			fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);
forrest@0: 			exit(4);
forrest@0: 		}
forrest@0: 		/* Get seconds of day variable */
forrest@0: 		wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
forrest@0: 		wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));
forrest@0: 		/* Get variable ID for requested variable */
forrest@0: 		wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
forrest@0: 		/* Retrieve time slice of desired variable */
forrest@0: 		for (ts = 0; ts < len; ts++) {
forrest@0: 			start[0] = ts;
forrest@0: 			count[0] = 1;
forrest@0: 			start[(ndims-2)] = 0;
forrest@0: 			count[(ndims-2)] = d1;
forrest@0: 			start[(ndims-1)] = 0;
forrest@0: 			count[(ndims-1)] = d2;
forrest@0: 			if (ndims == 4) {
forrest@0: 				start[1] = level;
forrest@0: 				count[1] = 1;
forrest@0: 			}
forrest@0: 			switch(nctype) {
forrest@0: 				case NC_FLOAT:
forrest@0: 					wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));
forrest@0: 					break;
forrest@0: 				case NC_DOUBLE:
forrest@0: 					wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));
forrest@0: 					break;
forrest@0: 				default:
forrest@0: 					fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 					exit(3);
forrest@0: 			}
forrest@0: 			nslice++;
forrest@0: 		}
forrest@0: 
forrest@0: 		/* Close input file */
forrest@0: 		wrap_nc(nc_close(input_ncid));
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: 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)
forrest@0: {
forrest@0: 	int i;
forrest@0: 
forrest@0: 	for (i = 0; i < nsamples; i++) {
forrest@0: 		switch(nctype) {
forrest@0: 			case NC_FLOAT:
forrest@0: 				if (stat_type == MEAN_TYPE)
forrest@0: 					add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0: 				else
forrest@0: 					add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0: 				break;
forrest@0: 			case NC_DOUBLE:
forrest@0: 				if (stat_type == MEAN_TYPE)
forrest@0: 					add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0: 				else
forrest@0: 					add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0: 				break;
forrest@0: 			default:
forrest@0: 				fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 				exit(3);
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	/* Divide sums by number of samples */
forrest@0: 	if (stat_type == MEAN_TYPE)
forrest@0: 		divide_means(d1, d2, means, cell_samples);
forrest@0: 	else
forrest@0: 		divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: float *double_to_float(size_t len, double *dvar)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	float *fvar;
forrest@0: 
forrest@0: 	if (!(fvar = (float *)malloc(len * sizeof(float)))) {
forrest@0: 		perror("double_to_float:fvar");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 
forrest@0: 	for (i = 0; i < len; i++)
forrest@0: 		fvar[i] = (float)dvar[i];
forrest@0: 
forrest@0: 	return fvar;
forrest@0: }
forrest@0: 
forrest@0: void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
forrest@0: 	int level, int ndims, double **var_hours)
forrest@0: {
forrest@0: 	/* Output dimensions should be one larger than input dimensions */
forrest@0: 	int i, hr;
forrest@0: 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 	float *fvar = NULL;
forrest@0: 
forrest@0: 	if (nctype == NC_FLOAT) {
forrest@0: 		if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
forrest@0: 			perror("write_var_hours:fvar");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0: 		start[0] = 0;
forrest@0: 		count[0] = 1; /* time */
forrest@0: 		start[1] = hr;
forrest@0: 		count[1] = 1; /* hour */
forrest@0: 		start[(ndims-2)] = 0;
forrest@0: 		count[(ndims-2)] = d1;
forrest@0: 		start[(ndims-1)] = 0;
forrest@0: 		count[(ndims-1)] = d2;
forrest@0: 		if (ndims == 5) {
forrest@0: 			start[2] = level;
forrest@0: 			count[2] = 1;
forrest@0: 		}
forrest@0: 		switch (nctype) {
forrest@0: 			case NC_FLOAT:
forrest@0: 				for (i = 0; i < (d1 * d2); i++)
forrest@0: 					fvar[i] = (float)var_hours[hr][i];
forrest@0: 				wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
forrest@0: 				break;
forrest@0: 			case NC_DOUBLE:
forrest@0: 				wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
forrest@0: 				break;
forrest@0: 			default:
forrest@0: 				fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0: 				exit(3);
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	if (nctype == NC_FLOAT)
forrest@0: 		free(fvar);
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
forrest@0: {
forrest@0: 	int j, k, nlevels, *mcsec;
forrest@0: 	size_t d1, d2, **cell_samples;
forrest@0: 	double **means, **stddevs;
forrest@0: 	struct var *in_vnode, *out_vnode;
forrest@0: 	void **vals;
forrest@0: 
forrest@0: 	if (input_var_head) {
forrest@0: 		in_vnode = input_var_head;
forrest@0: 		/* March through non-metadata variables to compute stats */
forrest@0: 		for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
forrest@0: 			if (!in_vnode->metadata) {
forrest@0: 				/* Make sure time is really the first dimension */
forrest@0: 				if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
forrest@0: 					fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
forrest@0: 					exit(3);
forrest@0: 				}
forrest@0: 				/* Figure out input dimensions */
forrest@0: 				if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
forrest@0: 					fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
forrest@0: 					exit(3);
forrest@0: 				}
forrest@0: 				/* Assume 2D output; find dimensions */
forrest@0: 				d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
forrest@0: 				d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
forrest@0: 				if (in_vnode->ndims == 3)
forrest@0: 					nlevels = 1;
forrest@0: 				else
forrest@0: 					nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
forrest@0: 				/* Allocate working space for means and stddevs */
forrest@0: 				alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
forrest@0: 				init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
forrest@0: 				/* Allocate and initize space for entire field across all files */
forrest@0: 				alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);
forrest@0: 				init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);
forrest@0: 				printf("Computing statistics for %s\n",
forrest@0: 					in_vnode->name);
forrest@0: 				/* Compute means and stddevs, then write them */
forrest@0: 				for (k = 0; k < nlevels; k++) {
forrest@0: 					/* Read the entire fields from all files */
forrest@0: 					read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);
forrest@0: 					/* Compute means */
forrest@0: 					compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
forrest@0: 					/* Find corresponding output variable on the mean output file */
forrest@0: 					out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
forrest@0: 					/* Write out the means for this variable */
forrest@0: 					write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
forrest@0: 					/* Zero out cell_samples so they can be used as a flag again for computing stddevs */
forrest@0: 					reset_cell_samples(d1, d2, cell_samples);
forrest@0: 					/* Compute stddevs using means */
forrest@0: 					compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
forrest@0: 					/* Find corresponding output variable on the stddev output file */
forrest@0: 					out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
forrest@0: 					/* Write out stddevs for this variable */
forrest@0: 					write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
forrest@0: 				}
forrest@0: 
forrest@0: 				/* Free all samples */
forrest@0: 				free_all_samples(nsamples, vals, mcsec);
forrest@0: 				/* Free working space for means and stddevs */
forrest@0: 				free_means_stddevs(means, stddevs, cell_samples);
forrest@0: 			}
forrest@0: 			in_vnode = in_vnode->next;
forrest@0: 		}
forrest@0: 	}
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void usage(char *program)
forrest@0: {
forrest@0: 	fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
forrest@0: 	fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
forrest@0: 	exit(4);
forrest@0: }
forrest@0: 
forrest@0: int main(int argc, char **argv)
forrest@0: {
forrest@0: 	int i, ic, nifnames;
forrest@0: 	size_t len, pos, nsamples;
forrest@0: 	char *mfname, *sfname, **ifnames, *flist;
forrest@0: 
forrest@0: 	ifnames = NULL;
forrest@0: 	mfname = sfname = flist = NULL;
forrest@0: 	nifnames = 0;
forrest@0: 
forrest@0: #ifdef DEBUG
forrest@0: 	printf("Number of metadata variables (nmvars) = %d\n", nmvars);
forrest@0: 	fflush(stdout);
forrest@0: #endif /* DEBUG */
forrest@0: 
forrest@0: 	/* check command-line flags and switches */
forrest@0: 	while ((ic = getopt(argc, argv, "m:s:")) != -1) {
forrest@0: 		switch(ic) {
forrest@0: 			case 'm':
forrest@0: 				if (!(mfname = strdup(optarg))) {
forrest@0: 					perror("mfname");
forrest@0: 					exit(2);
forrest@0: 				}
forrest@0: 				break;
forrest@0: 			case 's':
forrest@0: 				if (!(sfname = strdup(optarg))) {
forrest@0: 					perror("sfname");
forrest@0: 					exit(2);
forrest@0: 				}
forrest@0: 				break;
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	if (!mfname) {
forrest@0: 		fprintf(stderr, "Output file name for writing means required.\n");
forrest@0: 		usage(argv[0]);
forrest@0: 	}
forrest@0: 	if (!sfname) {
forrest@0: 		fprintf(stderr, "Output file name for writing standard deviations required.\n");
forrest@0: 		usage(argv[0]);
forrest@0: 	}
forrest@0: 
forrest@0: 	if (optind < argc) {
forrest@0: 		if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
forrest@0: 			perror("ifnames");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 		for (i = optind; i < argc; i++) {
forrest@0: 			if (!(ifnames[nifnames++] = strdup(argv[i]))) {
forrest@0: 				perror("ifnames[nifnames++]");
forrest@0: 				exit(2);
forrest@0: 			}
forrest@0: 		}
forrest@0: 	}
forrest@0: 	else {
forrest@0: 		fprintf(stderr, "At least one input file name is required.\n");
forrest@0: 		usage(argv[0]);
forrest@0: 	}
forrest@0: 
forrest@0: 
forrest@0: 	/*
forrest@0: 	 * Build list of input files to be included in the global history
forrest@0: 	 * attribute on the two outputs files.
forrest@0: 	 */
forrest@0: 	for (i = len = 0; i < nifnames; i++)
forrest@0: 		len += strlen(ifnames[i]);
forrest@0: 	len += nifnames + 1; /* space in front of every name + null terminator */
forrest@0: 	if (!(flist = (char *)malloc(len * sizeof(char)))) {
forrest@0: 		perror("flist");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
forrest@0: 		sprintf(&flist[pos], " %s", ifnames[i]);
forrest@0: #ifdef DEBUG
forrest@0: 	printf("flist=%s\n", flist);
forrest@0: 	fflush(stdout);
forrest@0: #endif /* DEBUG */
forrest@0: 
forrest@0: 	/* Open every file to count up the number of time samples. */
forrest@0: 	nsamples = count_samples(nifnames, ifnames);
forrest@0: #ifdef DEBUG
forrest@0: 	printf("Number of samples across all files = %ld\n", (long)nsamples);
forrest@0: 	fflush(stdout);
forrest@0: #endif /* DEBUG */
forrest@0: 
forrest@0: 	/*
forrest@0: 	 * Open the *last* input file and the two output files (for mean and
forrest@0: 	 * standard deviation).  Define dimensions, variables, and attributes
forrest@0: 	 * for the two output files based on the *last* input files.  The
forrest@0: 	 * last file is used because some metadata variables will contain
forrest@0: 	 * only values from the last time slice from the period over which
forrest@0: 	 * the statistics are computed.
forrest@0: 	 */
forrest@1: 	open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
forrest@1: 		nsamples);
forrest@0: 
forrest@0: 	compute_stats(nifnames, ifnames, nsamples);
forrest@0: 
forrest@0: 	/* Close the two output files */
forrest@0: 	wrap_nc(nc_close(mean_ncid));
forrest@0: 	wrap_nc(nc_close(stddev_ncid));
forrest@0: 
forrest@0: 	return 0;
forrest@0: }