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