#include <stdio.h>

#include <unistd.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

#include "netcdf.h"

/*

 * Loop through one month of hourly history tapes from the Community Land Model

 * and generate monthly statistics (means and standard deviations) of fields

 * by hour of day.

 */

#define	MEAN_TYPE	0

#define STDDEV_TYPE	1

#define NUM_TYPES	2

#define HOURS_PER_DAY	24

#define MIN_PER_HOUR	60

#define SEC_PER_MIN	60

#define SEC_PER_HOUR	(MIN_PER_HOUR * SEC_PER_MIN)

static char *metadata_vars[] = {

	"lon",

	"lat",

	"lonrof",

	"latrof",

	"time",

	"hour", /* new metadata variable to be added to output files */

	"levsoi",

	"levlak",

	"edgen",

	"edgee",

	"edges",

	"edgew",

	"longxy",

	"latixy",

	"area",

	"areaupsc",

	"topo",

	"topodnsc",

	"landfrac",

	"landmask",

	"pftmask",

	"indxupsc",

	"mcdate",

	"mcsec",

	"mdcur",

	"mscur",

	"nstep",

	"time_bounds",

	"date_written",

	"time_written"

};

struct text_att {

	char *name;

	char *value;

	struct text_att *next;

	struct text_att *prev;

};

struct dim {

	int dimid;

	char *name;

	size_t len;

	struct dim *next;

	struct dim *prev;

};

struct var {

	int ncvarid;

	char *name;

	nc_type nctype;

	int ndims;

	int *dimids;

	int natts;

	char metadata;

	char FillFlag;

	float FillValue;

	struct var *next;

	struct var *prev;

};

static char *time_name = "time";

static char *time_bounds_attname = "bounds";

static char *time_bounds_name = "time_bounds";

static char *climatology_bounds_attname = "climatology";

static char *climatology_bounds_name = "climatology_bounds";

static char *mcdate_name = "mcdate";

static char *mcsec_name = "mcsec";

static char *history_name = "history";

static char *nsamples_name = "nsamples";

static char *cell_method_name = "cell_method";

static char *cell_methods_prefix[] = { "time: mean within days time: mean over days", "time: mean within days times: standard_deviation over days" };

/* input stuff */

static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);

static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;

static struct text_att *input_gatt_head = NULL;

static struct dim *input_dim_head = NULL, **input_dim_idx;

static struct var *input_var_head = NULL;

/* translation stuff */

static int *idid2mdid, *idid2sdid; /* dimension IDs */

/* output stuff */

static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;

static struct dim *mean_dim_head = NULL;

static struct var *mean_var_head = NULL;

static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;

static struct dim *stddev_dim_head = NULL;

static struct var *stddev_var_head = NULL;

/* output values */

static double *times;

static double *tbounds;

static int *mcdate;

static int *mcsec;

char is_metadata(char *name)

{

	int i;

	for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);

	if (i < nmvars)

		return 1;

	else

		return 0;

}

#if 0

struct dim *dimlist_find_by_name(struct dim *head, char *name)

{

	int i;

	struct dim *p = NULL;

	if (head) {

		node = head;

		for (i = 0 ; i == 0 || node != head; i++) {

			if (!strcmp(node->name, name))

				p = node;

			node = node->next;

		}

		return p;

	}

	else

		return NULL;

}

#endif

struct var *varlist_find_by_name(struct var *head, char *name)

{

	int i;

	struct var *node;

	if (head) {

		node = head;

		for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)

			node = node->next;

		if (i && node == head)

			return NULL;

		else

			return node;

	}

	else

		return NULL;

}

void gattlist_add_node(struct text_att **headp, char *name, char *value)

{

	struct text_att *head, *node;

	head = *headp;

	if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {

		perror("gattlist_add_node");

		exit(2);

	}

	if (!(node->name = strdup(name))) {

		perror("gattlist_add_node");

		exit(2);

	}

	if (!(node->value = strdup(value))) {

		perror("gattlist_add_node");

		exit(2);

	}

	if (head == NULL) {

		/* first one */

		*headp = node;

		node->next = node;

		node->prev = node;

	}

	else {

		node->next = head;

		node->prev = head->prev;

		/* set this after setting node->prev from here */

		head->prev = node;

		/* set this after having set node->prev */

		node->prev->next = node;

	}

	return;

}

struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)

{

	struct dim *head, *node;

	head = *headp;

	if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {

		perror("dimlist_add_node");

		exit(2);

	}

	node->dimid = dimid;

	if (!(node->name = strdup(name))) {

		perror("dimlist_add_node");

		exit(2);

	}

	node->len = len;

	if (head == NULL) {

		/* first one */

		*headp = node;

		node->next = node;

		node->prev = node;

	}

	else {

		node->next = head;

		node->prev = head->prev;

		/* set this after setting node->prev from here */

		head->prev = node;

		/* set this after having set node->prev */

		node->prev->next = node;

	}

	return node;

}

void varlist_add_node(struct var **headp, int ncvarid, char *name,

	nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,

	float FillValue)

{

	int i;

	struct var *head, *node;

	head = *headp;

	if (!(node = (struct var *)malloc(sizeof(struct var)))) {

		perror("varlist_add_node");

		exit(3);

	}

	node->ncvarid = ncvarid;

	if (!(node->name = strdup(name))) {

		perror("varlist_add_node");

		exit(3);

	}

	node->nctype = nctype;

	node->ndims = ndims;

	if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {

		perror("varlist_add_node");

		exit(3);

	}

	for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];

	node->natts = natts;

	node->metadata = is_metadata(name);

	node->FillFlag = FillFlag;

	node->FillValue = FillValue;

	if (head == NULL) {

		/* first one */

		*headp = node;

		node->next = node;

		node->prev = node;

	}

	else {

		node->next = head;

		node->prev = head->prev;

		/* set this after setting node->prev from here */

		head->prev = node;

		/* set this after having set node->prev */

		node->prev->next = node;

	}

	return;

}

void varlist_print(struct var *headp)

{

	int i, j;

	struct var *node;

	printf("Beginning of Variable List\n");

	if (headp) {

		node = headp;

		for (j = 0; j == 0 || node != headp; j++) {

			printf("Variable %d (ptr=%ld):\n", j, (long)node);

			printf("\tncvarid = %d\n", node->ncvarid);

			printf("\tname = %s\n",  node->name);

			printf("\tnctype = %d\n", node->nctype);

			printf("\tndims = %d\n", node->ndims);

			printf("\tdimids =");

			for (i = 0; i < node->ndims; i++)

				printf(" %d", node->dimids[i]);

			printf("\n");

			printf("\tnatts = %d\n", node->natts);

			printf("\tmetadata = %d\n", (int)node->metadata);

			printf("\tnext = %ld\n", (long)node->next);

			printf("\tprev = %ld\n", (long)node->prev);

			node = node->next;

		}

	}

	else {

		printf("\tList undefined\n");

	}

	printf("End of Variable List\n");

	return;

}

void wrap_nc(int status)

{

	if (status != NC_NOERR) {

		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));

		exit(1);

	}

	return;

}

void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)

{

	int i;

	char name[NC_MAX_NAME+1];

	size_t len;

	struct dim **in_dim_idx;

	if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {

		perror("*in_dim_idxp");

		exit(3);

	}

	in_dim_idx = *in_dim_idxp;

	for (i = 0; i < ndims; i++) {

		wrap_nc(nc_inq_dim(ncid, i, name, &len));

		in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);

	}

	return;

}

void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)

{

	int i;

	char name[NC_MAX_NAME+1], *value;

	nc_type xtype;

	size_t len;

	for (i = 0; i < ngatts; i++) {

		wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));

		wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));

		if (xtype != NC_CHAR) {

			fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);

			exit(2);

		}

		if (!(value = (char *)malloc((len+1)*sizeof(char)))) {

			perror("get_gatts: value");

			exit(3);

		}

		wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));

		value[(len+1-1)] = '\0';

		gattlist_add_node(gatt_headp, name, value);

		free(value); /* because gattlist_add_node() duplicates it */

	}

	return;

}

void get_vars(int ncid, int nvars, struct var **var_headp)

{

	int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;

	size_t len;

	float FillValue;

	char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;

	nc_type xtype, atype;

	for (i = 0; i < nvars; i++) {

		wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,

			&natts));

		/* Check for _FillValue */

		FillFlag = 0;

		FillValue = 0.;

		if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)

			== NC_NOERR) {

			if (atype == NC_FLOAT && len) {

				wrap_nc(nc_get_att_float(ncid, i,

					fill_att_name, &FillValue));

				FillFlag = 1;

			}

		}

		varlist_add_node(var_headp, i, name, xtype, ndims, dimids,

			natts, FillFlag, FillValue);

	}

	return;

}

int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,

	size_t nsamples, int **in2outp, int out_ncid,

	struct dim **out_dim_headp, int *out_unlimdimidp)

{

	/*

	 * Define dimensions on new files and build translation tables between

	 * dimension IDs.

	 */

	int j, dimid, ndims, *in2out;

	size_t len;

	struct dim *dnode;

	ndims = 0;

	if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {

		perror("put_dimensions:in2outp");

		exit(3);

	}

	in2out = *in2outp;

	if (in_dim_head) {

		dnode = in_dim_head;

		for (j = 0; j == 0 || dnode != in_dim_head; j++) {

			if (dnode->dimid == in_unlimdimid)

				len = NC_UNLIMITED;

			else

				len = dnode->len;

			wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));

			dimlist_add_node(out_dim_headp, dimid, dnode->name, len);

			++ndims;

			in2out[dnode->dimid] = dimid;

			if (dnode->dimid == in_unlimdimid)

				*out_unlimdimidp = dimid;

			/* Just after the "time" dimension, add the new

			 * "nsamples" dimension. */

			if (!strcmp(dnode->name, time_name)) {

				wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));

				dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);

				++ndims;

			}

			dnode = dnode->next;

		}

	}

	else {

		fprintf(stderr, "WARNING: No dimensions defined!\n");

	}

	return ndims;

}

int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)

{

	/*

	 * Write out global attributes matching those of the input file.

	 * Change history attribute to the string passed in as an argument.

	 */

	int j, hflag = 0, ngatts;

	char *value;

	struct text_att *anode;

	ngatts = 0;

	if (in_gatt_head) {

		anode = in_gatt_head;

		for (j = 0; j == 0 || anode != in_gatt_head; j++) {

			if (!strcmp(anode->name, history_name)) {

				value = out_history;

				++hflag;

			}

			else

				value = anode->value;

			wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));

			++ngatts;

			anode = anode->next;

		}

		/* If no history attribute on input, add one to the output */

		if (!hflag) {

			wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));

			++ngatts;

		}

	}

	else {

		fprintf(stderr, "WARNING: No global attributes defined!\n");

	}

	return ngatts;

}

int translate_dimids(int in_ndims, int *in_dimids, int *in2out, int *out_dimids)

{

	/*

	 * Translate between input dimension IDs and those for the output file.

	 */

	int i, ndims;

	for (i = ndims = 0; i < in_ndims; i++, ndims++)

		out_dimids[ndims] = in2out[in_dimids[i]];

	return ndims;

}

int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,

	int in_varid, int out_ncid, int out_varid, char **cell_methods)

{

	/* 

	 * Copy the attributes of the input variable to those of the output

	 * variable. If the variable is not a metadata variable, ensure that

	 * the cell_method attribute is set correctly, depending on the output

	 * file type.

	 */

	int i, natts;

	char name[NC_MAX_NAME+1], cmflag = 0;

	for (i = natts = 0; i < in_natts; i++, natts++) {

		wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));

		if (!strcmp(name, cell_method_name) && !metadata) {

			wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));

			cmflag = 1;

		}

		else if (!strcmp(name, time_bounds_attname))

			/* Change time_bounds to climatology_bounds */

			wrap_nc(nc_put_att_text(out_ncid, out_varid, climatology_bounds_attname, strlen(climatology_bounds_name), climatology_bounds_name));

		else

			wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));

	}

	/*

	 * If no cell_method attribute was specified for a non-metadata

	 * variable on the input file, add the appropriate cell_method anyway

	 * on the output file.

	 */

	if (!cmflag && !metadata) {

		wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));

		natts++;

	}

	return natts;

}

int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,

	int *in2out, int out_ncid, struct var **out_var_headp,

	char **cell_methods, int stat_type)

{

	/*

	 * Define variables on output file and copy attributes from input file.

	 * Return number of variables defined.

	 */

	int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;

	struct var *vnode;

	nvars = 0;

	if (in_var_head) {

		vnode = in_var_head;

		/*

		 * March through input variables creating (mostly) the same

		 * variables on the output file.

		 */

		for (j = 0; j == 0 || vnode != in_var_head; j++) {

			ndims = translate_dimids(vnode->ndims, vnode->dimids, in2out, dimids);

			/* Create all normal data variables, but only the four

			 * time-based metadata variables: time,

			 * climatology_bounds, mcdate, and mcsec. */

			if (!strcmp(vnode->name, time_name)) {

				/* Magically promote "time" variable to

				 * double precision on output file */

				wrap_nc(nc_def_var(out_ncid, vnode->name, NC_DOUBLE, ndims, dimids, &varid));

				natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);

				varlist_add_node(out_var_headp, varid, vnode->name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);

				++nvars;

				/* Create "climatology_bounds" right after

				 * time and instead of "time_bounds", with no

				 * attributes */

				wrap_nc(nc_def_var(out_ncid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, &varid));

				natts = 0;

				varlist_add_node(out_var_headp, varid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);

				++nvars;

			}

			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))) {

				/* If it is not a metadata variable OR it is

				 * mcdate OR it is mcsec OR (it is a metadata

				 * variable that does not have time as its

				 * first dimension), then create it */

				wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));

				natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);

				varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);

				++nvars;

			}

			vnode = vnode->next;

		}

	}

	else {

		fprintf(stderr, "WARNING: No variables defined!\n");

	}

	return nvars;

}

void *alloc_var(nc_type nctype, size_t len)

{

	void *val;

	switch (nctype) {

		case NC_FLOAT:

			if (!(val = (float *)malloc((len) * sizeof(float)))) {

				perror("alloc_var: val");

				exit(3);

			}

			break;

		case NC_INT:

			if (!(val = (int *)malloc((len) * sizeof(int)))) {

				perror("alloc_var: val");

				exit(3);

			}

			break;

		case NC_DOUBLE:

			if (!(val = (double *)malloc((len) * sizeof(double)))) {

				perror("alloc_var: val");

				exit(3);

			}

			break;

		case NC_CHAR:

			if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {

				perror("alloc_var: val");

				exit(3);

			}

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return val;

}

void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)

{

	int i;

	size_t len = (size_t)1;

	void *val;

	/* Figure out the total size */

	for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;

	val = alloc_var(nctype, len);

	switch (nctype) {

		case NC_FLOAT:

			wrap_nc(nc_get_var_float(ncid, varid, val));

			break;

		case NC_INT:

			wrap_nc(nc_get_var_int(ncid, varid, val));

			break;

		case NC_DOUBLE:

			wrap_nc(nc_get_var_double(ncid, varid, val));

			break;

		case NC_CHAR:

			wrap_nc(nc_get_var_text(ncid, varid, val));

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return val;

}

void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)

{

	int i;

	size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

	void *val;

	/* Make sure time is really the first dimension */

	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {

		fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);

		exit(3);

	} 

	/*

	 * Figure out the total size for the slice (start at second dimension)

	 * and build start/count vectors.

	 */

	start[0] = tslice;

	count[0] = 1;

	for (i = 1; i < ndims; i++) {

		start[i] = 0;

		count[i] = (dim_idx[dimids[i]])->len;

		len *= (dim_idx[dimids[i]])->len;

	}

	val = alloc_var(nctype, len);

	switch (nctype) {

		case NC_FLOAT:

			wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));

			break;

		case NC_INT:

			wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));

			break;

		case NC_DOUBLE:

			wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));

			break;

		case NC_CHAR:

			wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return val;

}

void write_var(int ncid, int varid, nc_type nctype, void *val)

{

	switch (nctype) {

		case NC_FLOAT:

			wrap_nc(nc_put_var_float(ncid, varid, val));

			break;

		case NC_INT:

			wrap_nc(nc_put_var_int(ncid, varid, val));

			break;

		case NC_DOUBLE:

			wrap_nc(nc_put_var_double(ncid, varid, val));

			break;

		case NC_CHAR:

			wrap_nc(nc_put_var_text(ncid, varid, val));

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return;

}

void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)

{

	int i;

	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

	/* Make sure time is really the first dimension */

	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {

		fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);

		exit(3);

	}

	/* Build start/count vectors */

	start[0] = tslice;

	count[0] = 1;

	for (i = 1; i < ndims; i++) {

		start[i] = 0;

		count[i] = (dim_idx[dimids[i]])->len;

	}

	switch (nctype) {

		case NC_FLOAT:

			wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));

			break;

		case NC_INT:

			wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));

			break;

		case NC_DOUBLE:

			wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));

			break;

		case NC_CHAR:

			wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return;

}

void write_time_var(int ncid, int varid, nc_type nctype, void *time_var)

{

	size_t start[1], count[1];

	start[0] = 0;

	count[0] = 24;

	switch (nctype) {

		case NC_FLOAT:

			wrap_nc(nc_put_vara_float(ncid, varid, start, count, time_var));

			break;

		case NC_INT:

			wrap_nc(nc_put_vara_int(ncid, varid, start, count, time_var));

			break;

		case NC_DOUBLE:

			wrap_nc(nc_put_vara_double(ncid, varid, start, count, time_var));

			break;

		case NC_CHAR:

			wrap_nc(nc_put_vara_text(ncid, varid, start, count, time_var));

			break;

		default:

			fprintf(stderr, "netCDF external data type %d not supported\n", nctype);

			exit(3);

	}

	return;

}

void write_tbounds(int ncid, int varid, double *tbounds)

{

	size_t start[2], count[2];

	start[0] = start[1] = 0;

	count[0] = 24; count[1] = 2;

	wrap_nc(nc_put_vara_double(ncid, varid, start, count, tbounds));

	return;

}

void copy_metadata(int in_ncid, struct var *in_var_head,

	struct dim **in_dim_idx, int out_ncid, struct var *out_var_head,

	double *times, double *tbounds, int *mcdate, int *mcsec)

{

	int i;

	struct var *in_vnode, *out_vnode;

	void *val;

	if (!in_var_head) return;

	in_vnode = in_var_head;

	/*

	 * March through input variables to stuff metadata values into

	 * the new output files. NOTE: The only time-based metadata

	 * variables, explicitly handled, are: time, climatology_bounds,

	 * mcdate, and mcsec.  These are passed in.

	 */

	for (i = 0; i == 0 || in_vnode != in_var_head; i++) {

		if (in_vnode->metadata) {

			/* Find the corresponding output variable */

			if (!strcmp(in_vnode->name, time_bounds_name))

				out_vnode = varlist_find_by_name(out_var_head, climatology_bounds_name);

			else

				out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);

			if (out_vnode) {

				/* Read values from input and write to output */

				if (strcmp((in_dim_idx[in_vnode->dimids[0]])->name, time_name)) {

					/* "time" is not the first dimension,

					 * so just copy the variable values */

#ifdef DEBUG

					printf("Copying metadata variable %s\n",

						in_vnode->name);

#endif /* DEBUG */

					val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);

					write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);

					free(val);

				}

				else {

					/* "time" is the first dimension, so

					 * explicitly handle the "time",

					 * "climatology_bounds", "mcdate", and

					 * "mcsec" while dropping all others */

					if (!strcmp(out_vnode->name, time_name)) {

						/* Force time stamps to be

						 * those computed previously,

						 * and remember these are now

						 * double precision. */

#ifdef DEBUG

						printf("Writing metadata variable %s\n", out_vnode->name);

#endif /* DEBUG */

						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, times);

					}

					if (!strcmp(out_vnode->name, climatology_bounds_name)) {

						/* time_bounds are now

						 * climatology_bounds and have

						 * pre-computed values */

#ifdef DEBUG

						printf("Writing metadata variable %s\n", out_vnode->name);

#endif /* DEBUG */

						write_tbounds(out_ncid, out_vnode->ncvarid, tbounds);

					}

					else if (!strcmp(out_vnode->name, mcdate_name)) {

						/* mcdate is pre-computed */

#ifdef DEBUG

						printf("Writing metadata variable %s\n", out_vnode->name);

#endif /* DEBUG */

						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcdate);

					}

					else if (!strcmp(out_vnode->name, mcsec_name)) {

						/* mcsec is pre-computed */

#ifdef DEBUG

						printf("Writing metadata variable %s\n", out_vnode->name);

#endif /* DEBUG */

						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcsec);

					}

				}

			}

			else {

#ifdef DEBUG

				printf("Dropping metadata variable %s\n", in_vnode->name);

#endif /* DEBUG */

			}

		}

		else {

#ifdef DEBUG

			printf("Skipping non-metadata variable %s\n",

				in_vnode->name);

#endif /* DEBUG */

		}

#ifdef DEBUG

		printf ("Done\n");

#endif /* DEBUG */

		in_vnode = in_vnode->next;

	}

	return;

}

void get_time_bounds(char *first_fname, char *last_fname, double *times, double *tbounds, int *mcdate, int *mcsec)

{

	int i, time_dimid, time_bounds_varid, mcdate_varid, mcsec_varid;

	size_t time_len, time_bounds_start[2], time_bounds_count[2],

		mc_start[1], mc_count[1];

	double bnd1[24], bnd2[24], day1, day2;

	/* Open first input file */

	wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));

	/* Get dimension ID of the time dimension */

	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));

	/* Get length of time dimension */

	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));

	/* Get variable ID of time_bounds variable */

	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));

	/* Set start/count for reading the left values of time_bounds from the

	 * first 24 timeslices of the first file */

	time_bounds_start[0] = time_bounds_start[1] = 0;

	time_bounds_count[0] = 24; time_bounds_count[1] = 1;

	/* Read the values */

	wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,

		time_bounds_start, time_bounds_count, bnd1));

	/* If the first and last file are not the same, close the first one and

	 * open the second one */

	if (strcmp(first_fname, last_fname)) {

		/* Close the first input file */

		wrap_nc(nc_close(input_ncid));

		/* Open the last input file */

		wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));

		/* Get dimension ID of the time dimension */

		wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));

		/* Get length of time dimension */

		wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));

		/* Get variable ID of time_bounds variable */

		wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,

			&time_bounds_varid));

	}

	/* Set start/count for reading the right values of time_bounds from the

	 * last 24 timeslices of the last file */

	time_bounds_start[0] = time_len - 24; time_bounds_start[1] = 1;

	time_bounds_count[0] = 24; time_bounds_count[1] = 1;

	/* Read the value */

	wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,

		time_bounds_start, time_bounds_count, bnd2));

	/* Read the last 24 mcdate and mcsec values */

	wrap_nc(nc_inq_varid(input_ncid, mcdate_name, &mcdate_varid));

	wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &mcsec_varid));

	mc_start[0] = time_len - 24;

	mc_count[0] = 24;

	wrap_nc(nc_get_vara_int(input_ncid, mcdate_varid, mc_start, mc_count,

		mcdate));

	wrap_nc(nc_get_vara_int(input_ncid, mcsec_varid, mc_start, mc_count,

		mcsec));

	/* Close the last input file */

	wrap_nc(nc_close(input_ncid));

	/* Store time bounds and compute new time stamps as midpoints */

	for (i = 0; i < 24; i++) {

		tbounds[(i*2)] = bnd1[i];

		tbounds[(i*2+1)] = bnd2[i];

		/* Pick one (integer) day near the midpoint of the bounds and

		 * add in the fraction from the right bounds, so for January

		 * the first hour-of-day value would be 15.0 (plus some year

		 * offset) */

		day1 = floor(bnd1[i]);

		day2 = floor(bnd2[i]);

		times[i] = floor((day1 + day2) / 2.0) + (bnd2[i]-day2);