c-lamp_postproc: h1_summary2.c@2ce4ee911439 (annotated)

forrest@0	1	#include <stdio.h>
forrest@0	2	#include <unistd.h>
forrest@0	3	#include <stdlib.h>
forrest@0	4	#include <string.h>
forrest@0	5	#include <math.h>
forrest@0	6	#include "netcdf.h"
forrest@0	7
forrest@0	8	/*
forrest@0	9	* Loop through one month of hourly history tapes from the Community Land Model
forrest@0	10	* and generate monthly statistics (means and standard deviations) of fields
forrest@0	11	* by hour of day.
forrest@0	12	*/
forrest@0	13
forrest@0	14	#define MEAN_TYPE 0
forrest@0	15	#define STDDEV_TYPE 1
forrest@0	16
forrest@0	17	#define HOURS_PER_DAY 24
forrest@0	18	#define MIN_PER_HOUR 60
forrest@0	19	#define SEC_PER_MIN 60
forrest@0	20	#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
forrest@0	21
forrest@0	22	static char *metadata_vars[] = {
forrest@0	23	"lon",
forrest@0	24	"lat",
forrest@0	25	"lonrof",
forrest@0	26	"latrof",
forrest@0	27	"time",
forrest@0	28	"hour", /* new metadata variable to be added to output files */
forrest@0	29	"levsoi",
forrest@0	30	"levlak",
forrest@0	31	"edgen",
forrest@0	32	"edgee",
forrest@0	33	"edges",
forrest@0	34	"edgew",
forrest@0	35	"longxy",
forrest@0	36	"latixy",
forrest@0	37	"area",
forrest@0	38	"areaupsc",
forrest@0	39	"topo",
forrest@0	40	"topodnsc",
forrest@0	41	"landfrac",
forrest@0	42	"landmask",
forrest@0	43	"pftmask",
forrest@0	44	"indxupsc",
forrest@0	45	"mcdate",
forrest@0	46	"mcsec",
forrest@0	47	"mdcur",
forrest@0	48	"mscur",
forrest@0	49	"nstep",
forrest@0	50	"time_bounds",
forrest@0	51	"date_written",
forrest@0	52	"time_written"
forrest@0	53	};
forrest@0	54
forrest@0	55	struct text_att {
forrest@0	56	char *name;
forrest@0	57	char *value;
forrest@0	58	struct text_att *next;
forrest@0	59	struct text_att *prev;
forrest@0	60	};
forrest@0	61
forrest@0	62	struct dim {
forrest@0	63	int dimid;
forrest@0	64	char *name;
forrest@0	65	size_t len;
forrest@0	66	struct dim *next;
forrest@0	67	struct dim *prev;
forrest@0	68	};
forrest@0	69
forrest@0	70	struct var {
forrest@0	71	int ncvarid;
forrest@0	72	char *name;
forrest@0	73	nc_type nctype;
forrest@0	74	int ndims;
forrest@0	75	int *dimids;
forrest@0	76	int natts;
forrest@0	77	char metadata;
forrest@0	78	char FillFlag;
forrest@0	79	float FillValue;
forrest@0	80	struct var *next;
forrest@0	81	struct var *prev;
forrest@0	82	};
forrest@0	83
forrest@0	84	static char *time_name = "time";
forrest@1	85	static char *time_bounds_name = "time_bounds";
forrest@0	86	static char *mcsec_name = "mcsec";
forrest@0	87	static char *history_name = "history";
forrest@0	88	static char *nsamples_name = "nsamples";
forrest@0	89	static char hour_name = "hour", hour_long_name = "hour of day",
forrest@0	90	*hour_units = "hour";
forrest@0	91	static char *cell_method_name = "cell_method";
forrest@0	92	/* input stuff */
forrest@0	93	static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
forrest@0	94	static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
forrest@0	95	static struct text_att *input_gatt_head = NULL;
forrest@0	96	static struct dim input_dim_head = NULL, *input_dim_idx;
forrest@0	97	static struct var *input_var_head = NULL;
forrest@0	98	/* translation stuff */
forrest@0	99	static int idid2mdid, idid2sdid; /* dimension IDs */
forrest@0	100	/* output stuff */
forrest@0	101	static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
forrest@0	102	static int mean_hour_dimid; /* special notes */
forrest@0	103	static struct dim *mean_dim_head = NULL;
forrest@0	104	static struct var *mean_var_head = NULL;
forrest@0	105	static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
forrest@0	106	static int stddev_hour_dimid; /* special notes */
forrest@0	107	static struct dim *stddev_dim_head = NULL;
forrest@0	108	static struct var *stddev_var_head = NULL;
forrest@0	109
forrest@0	110	char is_metadata(char *name)
forrest@0	111	{
forrest@0	112	int i;
forrest@0	113
forrest@0	114	for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
forrest@0	115
forrest@0	116	if (i < nmvars)
forrest@0	117	return 1;
forrest@0	118	else
forrest@0	119	return 0;
forrest@0	120	}
forrest@0	121
forrest@0	122	#if 0
forrest@0	123	struct dim dimlist_find_by_name(struct dim head, char *name)
forrest@0	124	{
forrest@0	125	int i;
forrest@0	126	struct dim *p = NULL;
forrest@0	127
forrest@0	128	if (head) {
forrest@0	129	node = head;
forrest@0	130	for (i = 0 ; i == 0 \|\| node != head; i++) {
forrest@0	131	if (!strcmp(node->name, name))
forrest@0	132	p = node;
forrest@0	133	node = node->next;
forrest@0	134	}
forrest@0	135	return p;
forrest@0	136	}
forrest@0	137	else
forrest@0	138	return NULL;
forrest@0	139	}
forrest@0	140	#endif
forrest@0	141
forrest@0	142	struct var varlist_find_by_name(struct var head, char *name)
forrest@0	143	{
forrest@0	144	int i;
forrest@0	145	struct var *node;
forrest@0	146
forrest@0	147	if (head) {
forrest@0	148	node = head;
forrest@0	149	for (i = 0 ; (i == 0 \|\| node != head) && strcmp(node->name, name); i++)
forrest@0	150	node = node->next;
forrest@0	151	if (i && node == head)
forrest@0	152	return NULL;
forrest@0	153	else
forrest@0	154	return node;
forrest@0	155	}
forrest@0	156	else
forrest@0	157	return NULL;
forrest@0	158	}
forrest@0	159
forrest@0	160	void gattlist_add_node(struct text_att *headp, char name, char *value)
forrest@0	161	{
forrest@0	162	struct text_att head, node;
forrest@0	163
forrest@0	164	head = *headp;
forrest@0	165
forrest@0	166	if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
forrest@0	167	perror("gattlist_add_node");
forrest@0	168	exit(2);
forrest@0	169	}
forrest@0	170
forrest@0	171	if (!(node->name = strdup(name))) {
forrest@0	172	perror("gattlist_add_node");
forrest@0	173	exit(2);
forrest@0	174	}
forrest@0	175	if (!(node->value = strdup(value))) {
forrest@0	176	perror("gattlist_add_node");
forrest@0	177	exit(2);
forrest@0	178	}
forrest@0	179
forrest@0	180	if (head == NULL) {
forrest@0	181	/* first one */
forrest@0	182	*headp = node;
forrest@0	183	node->next = node;
forrest@0	184	node->prev = node;
forrest@0	185	}
forrest@0	186	else {
forrest@0	187	node->next = head;
forrest@0	188	node->prev = head->prev;
forrest@0	189	/* set this after setting node->prev from here */
forrest@0	190	head->prev = node;
forrest@0	191	/* set this after having set node->prev */
forrest@0	192	node->prev->next = node;
forrest@0	193	}
forrest@0	194
forrest@0	195	return;
forrest@0	196	}
forrest@0	197
forrest@0	198	struct dim dimlist_add_node(struct dim headp, int dimid, char name, size_t len)
forrest@0	199	{
forrest@0	200	struct dim head, node;
forrest@0	201
forrest@0	202	head = *headp;
forrest@0	203
forrest@0	204	if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
forrest@0	205	perror("dimlist_add_node");
forrest@0	206	exit(2);
forrest@0	207	}
forrest@0	208
forrest@0	209	node->dimid = dimid;
forrest@0	210	if (!(node->name = strdup(name))) {
forrest@0	211	perror("dimlist_add_node");
forrest@0	212	exit(2);
forrest@0	213	}
forrest@0	214	node->len = len;
forrest@0	215
forrest@0	216	if (head == NULL) {
forrest@0	217	/* first one */
forrest@0	218	*headp = node;
forrest@0	219	node->next = node;
forrest@0	220	node->prev = node;
forrest@0	221	}
forrest@0	222	else {
forrest@0	223	node->next = head;
forrest@0	224	node->prev = head->prev;
forrest@0	225	/* set this after setting node->prev from here */
forrest@0	226	head->prev = node;
forrest@0	227	/* set this after having set node->prev */
forrest@0	228	node->prev->next = node;
forrest@0	229	}
forrest@0	230
forrest@0	231	return node;
forrest@0	232	}
forrest@0	233
forrest@0	234	void varlist_add_node(struct var *headp, int ncvarid, char name,
forrest@0	235	nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
forrest@0	236	float FillValue)
forrest@0	237	{
forrest@0	238	int i;
forrest@0	239	struct var head, node;
forrest@0	240
forrest@0	241	head = *headp;
forrest@0	242
forrest@0	243	if (!(node = (struct var *)malloc(sizeof(struct var)))) {
forrest@0	244	perror("varlist_add_node");
forrest@0	245	exit(3);
forrest@0	246	}
forrest@0	247
forrest@0	248	node->ncvarid = ncvarid;
forrest@0	249	if (!(node->name = strdup(name))) {
forrest@0	250	perror("varlist_add_node");
forrest@0	251	exit(3);
forrest@0	252	}
forrest@0	253	node->nctype = nctype;
forrest@0	254	node->ndims = ndims;
forrest@0	255	if (!(node->dimids = (int )malloc(ndims sizeof(int)))) {
forrest@0	256	perror("varlist_add_node");
forrest@0	257	exit(3);
forrest@0	258	}
forrest@0	259	for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
forrest@0	260	node->natts = natts;
forrest@0	261	node->metadata = is_metadata(name);
forrest@0	262	node->FillFlag = FillFlag;
forrest@0	263	node->FillValue = FillValue;
forrest@0	264
forrest@0	265	if (head == NULL) {
forrest@0	266	/* first one */
forrest@0	267	*headp = node;
forrest@0	268	node->next = node;
forrest@0	269	node->prev = node;
forrest@0	270	}
forrest@0	271	else {
forrest@0	272	node->next = head;
forrest@0	273	node->prev = head->prev;
forrest@0	274	/* set this after setting node->prev from here */
forrest@0	275	head->prev = node;
forrest@0	276	/* set this after having set node->prev */
forrest@0	277	node->prev->next = node;
forrest@0	278	}
forrest@0	279
forrest@0	280	return;
forrest@0	281	}
forrest@0	282
forrest@0	283	void varlist_print(struct var *headp)
forrest@0	284	{
forrest@0	285	int i, j;
forrest@0	286	struct var *node;
forrest@0	287
forrest@0	288	printf("Beginning of Variable List\n");
forrest@0	289
forrest@0	290	if (headp) {
forrest@0	291	node = headp;
forrest@0	292	for (j = 0; j == 0 \|\| node != headp; j++) {
forrest@0	293	printf("Variable %d (ptr=%ld):\n", j, (long)node);
forrest@0	294	printf("\tncvarid = %d\n", node->ncvarid);
forrest@0	295	printf("\tname = %s\n", node->name);
forrest@0	296	printf("\tnctype = %d\n", node->nctype);
forrest@0	297	printf("\tndims = %d\n", node->ndims);
forrest@0	298	printf("\tdimids =");
forrest@0	299	for (i = 0; i < node->ndims; i++)
forrest@0	300	printf(" %d", node->dimids[i]);
forrest@0	301	printf("\n");
forrest@0	302	printf("\tnatts = %d\n", node->natts);
forrest@0	303	printf("\tmetadata = %d\n", (int)node->metadata);
forrest@0	304	printf("\tnext = %ld\n", (long)node->next);
forrest@0	305	printf("\tprev = %ld\n", (long)node->prev);
forrest@0	306	node = node->next;
forrest@0	307	}
forrest@0	308	}
forrest@0	309	else {
forrest@0	310	printf("\tList undefined\n");
forrest@0	311	}
forrest@0	312
forrest@0	313	printf("End of Variable List\n");
forrest@0	314
forrest@0	315	return;
forrest@0	316	}
forrest@0	317
forrest@0	318	void wrap_nc(int status)
forrest@0	319	{
forrest@0	320	if (status != NC_NOERR) {
forrest@0	321	fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
forrest@0	322	exit(1);
forrest@0	323	}
forrest@0	324
forrest@0	325	return;
forrest@0	326	}
forrest@0	327
forrest@0	328	void get_dimensions(int ncid, int ndims, struct dim dim_headp, struct dim *in_dim_idxp)
forrest@0	329	{
forrest@0	330	int i;
forrest@0	331	char name[NC_MAX_NAME+1];
forrest@0	332	size_t len;
forrest@0	333	struct dim **in_dim_idx;
forrest@0	334
forrest@0	335	if (!(in_dim_idxp = (struct dim )malloc(ndims sizeof(struct dim *)))) {
forrest@0	336	perror("*in_dim_idxp");
forrest@0	337	exit(3);
forrest@0	338	}
forrest@0	339	in_dim_idx = *in_dim_idxp;
forrest@0	340
forrest@0	341	for (i = 0; i < ndims; i++) {
forrest@0	342	wrap_nc(nc_inq_dim(ncid, i, name, &len));
forrest@0	343	in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
forrest@0	344	}
forrest@0	345
forrest@0	346	return;
forrest@0	347	}
forrest@0	348	void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
forrest@0	349	{
forrest@0	350	int i;
forrest@0	351	char name[NC_MAX_NAME+1], *value;
forrest@0	352	nc_type xtype;
forrest@0	353	size_t len;
forrest@0	354
forrest@0	355	for (i = 0; i < ngatts; i++) {
forrest@0	356	wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
forrest@0	357	wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
forrest@0	358	if (xtype != NC_CHAR) {
forrest@0	359	fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
forrest@0	360	exit(2);
forrest@0	361	}
forrest@0	362	if (!(value = (char )malloc((len+1)sizeof(char)))) {
forrest@0	363	perror("get_gatts: value");
forrest@0	364	exit(3);
forrest@0	365	}
forrest@0	366	wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
forrest@0	367	value[(len+1-1)] = '\0';
forrest@0	368	gattlist_add_node(gatt_headp, name, value);
forrest@0	369	free(value); /* because gattlist_add_node() duplicates it */
forrest@0	370	}
forrest@0	371
forrest@0	372	return;
forrest@0	373	}
forrest@0	374
forrest@0	375	void get_vars(int ncid, int nvars, struct var **var_headp)
forrest@0	376	{
forrest@0	377	int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
forrest@0	378	size_t len;
forrest@0	379	float FillValue;
forrest@0	380	char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
forrest@0	381	nc_type xtype, atype;
forrest@0	382
forrest@0	383	for (i = 0; i < nvars; i++) {
forrest@0	384	wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
forrest@0	385	&natts));
forrest@0	386	/* Check for _FillValue */
forrest@0	387	FillFlag = 0;
forrest@0	388	FillValue = 0.;
forrest@0	389	if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
forrest@0	390	== NC_NOERR) {
forrest@0	391	if (atype == NC_FLOAT && len) {
forrest@0	392	wrap_nc(nc_get_att_float(ncid, i,
forrest@0	393	fill_att_name, &FillValue));
forrest@0	394	FillFlag = 1;
forrest@0	395	}
forrest@0	396	}
forrest@0	397	varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
forrest@0	398	natts, FillFlag, FillValue);
forrest@0	399	}
forrest@0	400
forrest@0	401	return;
forrest@0	402	}
forrest@0	403
forrest@0	404	int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
forrest@0	405	size_t nsamples, int **in2outp, int out_ncid,
forrest@0	406	struct dim *out_dim_headp, int out_unlimdimidp, int *out_hour_dimidp)
forrest@0	407	{
forrest@0	408	/*
forrest@0	409	* Define dimensions on new files and build translation tables between
forrest@0	410	* dimension IDs.
forrest@0	411	*/
forrest@0	412	int j, dimid, ndims, *in2out;
forrest@0	413	size_t len;
forrest@0	414	struct dim *dnode;
forrest@0	415
forrest@0	416	ndims = 0;
forrest@0	417
forrest@0	418	if (!(in2outp = (int )malloc((in_ndims+1)*sizeof(int)))) {
forrest@0	419	perror("put_dimensions: in2outp");
forrest@0	420	exit(3);
forrest@0	421	}
forrest@0	422	in2out = *in2outp;
forrest@0	423
forrest@0	424	if (in_dim_head) {
forrest@0	425	dnode = in_dim_head;
forrest@0	426	for (j = 0; j == 0 \|\| dnode != in_dim_head; j++) {
forrest@0	427	if (dnode->dimid == in_unlimdimid)
forrest@0	428	len = NC_UNLIMITED;
forrest@0	429	else
forrest@0	430	len = dnode->len;
forrest@0	431	wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
forrest@0	432	dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
forrest@0	433	++ndims;
forrest@0	434	in2out[dnode->dimid] = dimid;
forrest@0	435	if (dnode->dimid == in_unlimdimid)
forrest@0	436	*out_unlimdimidp = dimid;
forrest@0	437	/*
forrest@0	438	* Just after the "time" dimension, add the new
forrest@0	439	* "nsamples" and "hour" dimensions.
forrest@0	440	*/
forrest@0	441	if (!strcmp(dnode->name, time_name)) {
forrest@0	442	wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
forrest@0	443	dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
forrest@0	444	++ndims;
forrest@0	445
forrest@0	446	wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
forrest@0	447	dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
forrest@0	448	++ndims;
forrest@0	449	/* Track hour dimid for out files */
forrest@0	450	*out_hour_dimidp = dimid;
forrest@0	451	}
forrest@0	452
forrest@0	453	dnode = dnode->next;
forrest@0	454	}
forrest@0	455	}
forrest@0	456	else {
forrest@0	457	fprintf(stderr, "WARNING: No dimensions defined!\n");
forrest@0	458	}
forrest@0	459
forrest@0	460	return ndims;
forrest@0	461	}
forrest@0	462
forrest@0	463	int put_gatts(struct text_att in_gatt_head, int out_ncid, char out_history)
forrest@0	464	{
forrest@0	465	/*
forrest@0	466	* Write out global attributes matching those of the input file.
forrest@0	467	* Change history attribute to the string passed in as an argument.
forrest@0	468	*/
forrest@0	469	int j, hflag = 0, ngatts;
forrest@0	470	char *value;
forrest@0	471	struct text_att *anode;
forrest@0	472
forrest@0	473	ngatts = 0;
forrest@0	474
forrest@0	475	if (in_gatt_head) {
forrest@0	476	anode = in_gatt_head;
forrest@0	477	for (j = 0; j == 0 \|\| anode != in_gatt_head; j++) {
forrest@0	478	if (!strcmp(anode->name, history_name)) {
forrest@0	479	value = out_history;
forrest@0	480	++hflag;
forrest@0	481	}
forrest@0	482	else
forrest@0	483	value = anode->value;
forrest@0	484	wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
forrest@0	485	++ngatts;
forrest@0	486	anode = anode->next;
forrest@0	487	}
forrest@0	488	/* If no history attribute on input, add one to the output */
forrest@0	489	if (!hflag) {
forrest@0	490	wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
forrest@0	491	++ngatts;
forrest@0	492	}
forrest@0	493	}
forrest@0	494	else {
forrest@0	495	fprintf(stderr, "WARNING: No global attributes defined!\n");
forrest@0	496	}
forrest@0	497
forrest@0	498	return ngatts;
forrest@0	499	}
forrest@0	500
forrest@0	501	int translate_dimids(struct dim *in_dim_idx, char metadata, int in_ndims, int in_dimids, int in2out, int out_dimids, int hour_dimid)
forrest@0	502	{
forrest@0	503	/*
forrest@0	504	* Translate between input dimension IDs and those for the output file.
forrest@0	505	* For normal time-based variables, add a new dimension for hour of
forrest@0	506	* day. For metadata variables, do not add new dimensions, even if
forrest@0	507	* it is time-based. Return (possibly new) number of dimensions.
forrest@0	508	*/
forrest@0	509	int i, ndims;
forrest@0	510
forrest@0	511	for (i = ndims = 0; i < in_ndims; i++, ndims++) {
forrest@0	512	out_dimids[ndims] = in2out[in_dimids[i]];
forrest@0	513	if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
forrest@0	514	&& !metadata) {
forrest@0	515	ndims++;
forrest@0	516	out_dimids[ndims] = hour_dimid;
forrest@0	517	}
forrest@0	518	}
forrest@0	519
forrest@0	520	return ndims;
forrest@0	521	}
forrest@0	522
forrest@0	523	int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
forrest@0	524	int in_varid, int out_ncid, int out_varid)
forrest@0	525	{
forrest@0	526	/*
forrest@0	527	* Copy the attributes of the input variable to those of the output
forrest@0	528	* variable. If the variable is not a metadata variable, ensure that
forrest@0	529	* the cell_method attribute is set either to "time: mean" or
forrest@0	530	* "time: standard_deviation", depending on the output file type.
forrest@0	531	*/
forrest@0	532
forrest@0	533	int i, natts;
forrest@0	534	char name[NC_MAX_NAME+1], cmflag = 0;
forrest@0	535	char *cell_methods[] = { "time: mean", "time: standard_deviation" };
forrest@0	536
forrest@0	537	for (i = natts = 0; i < in_natts; i++, natts++) {
forrest@0	538	wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
forrest@0	539	if (!strcmp(name, cell_method_name) && !metadata) {
forrest@0	540	wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
forrest@0	541	cmflag = 1;
forrest@0	542	}
forrest@0	543	else
forrest@0	544	wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
forrest@0	545	}
forrest@0	546	/*
forrest@0	547	* If no cell_method attribute was specified for a non-metadata
forrest@0	548	* variable on the input file, add the appropriate cell_method anyway
forrest@0	549	* on the output file.
forrest@0	550	*/
forrest@0	551	if (!cmflag && !metadata) {
forrest@0	552	wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
forrest@0	553	natts++;
forrest@0	554	}
forrest@0	555
forrest@0	556	return natts;
forrest@0	557	}
forrest@0	558
forrest@0	559	int define_vars(int in_ncid, struct dim *in_dim_idx, struct var in_var_head,
forrest@0	560	int in2out, int out_ncid, int hour_dimid, struct var *out_var_headp,
forrest@0	561	int stat_type)
forrest@0	562	{
forrest@0	563	/*
forrest@0	564	* Define variables on output file and copy attributes from input file.
forrest@0	565	* Return number of variables defined.
forrest@0	566	*/
forrest@0	567	int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
forrest@0	568	struct var *vnode;
forrest@0	569
forrest@0	570	nvars = 0;
forrest@0	571
forrest@0	572	if (in_var_head) {
forrest@0	573	vnode = in_var_head;
forrest@0	574	/*
forrest@0	575	* March through input variables creating (mostly) the same
forrest@0	576	* variables on the output file.
forrest@0	577	*/
forrest@0	578	for (j = 0; j == 0 \|\| vnode != in_var_head; j++) {
forrest@0	579	ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
forrest@0	580	wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
forrest@0	581	natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
forrest@0	582	varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
forrest@0	583	++nvars;
forrest@0	584	/*
forrest@0	585	* Just after the "time" variable, add the new "hour"
forrest@0	586	* variable for hour of day.
forrest@0	587	*/
forrest@0	588	if (!strcmp(vnode->name, time_name)) {
forrest@0	589	ndims = 1;
forrest@0	590	dimids[0] = hour_dimid;
forrest@0	591	wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
forrest@0	592	wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
forrest@0	593	wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
forrest@0	594	varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
forrest@0	595	++nvars;
forrest@0	596	}
forrest@0	597
forrest@0	598	vnode = vnode->next;
forrest@0	599	}
forrest@0	600	}
forrest@0	601	else {
forrest@0	602	fprintf(stderr, "WARNING: No variables defined!\n");
forrest@0	603	}
forrest@0	604
forrest@0	605	return nvars;
forrest@0	606	}
forrest@0	607
forrest@0	608	void *alloc_var(nc_type nctype, size_t len)
forrest@0	609	{
forrest@0	610	void *val;
forrest@0	611
forrest@0	612	switch (nctype) {
forrest@0	613	case NC_FLOAT:
forrest@0	614	if (!(val = (float )malloc((len) sizeof(float)))) {
forrest@0	615	perror("alloc_var: val");
forrest@0	616	exit(3);
forrest@0	617	}
forrest@0	618	break;
forrest@0	619	case NC_INT:
forrest@0	620	if (!(val = (int )malloc((len) sizeof(int)))) {
forrest@0	621	perror("alloc_var: val");
forrest@0	622	exit(3);
forrest@0	623	}
forrest@0	624	break;
forrest@0	625	case NC_DOUBLE:
forrest@0	626	if (!(val = (double )malloc((len) sizeof(double)))) {
forrest@0	627	perror("alloc_var: val");
forrest@0	628	exit(3);
forrest@0	629	}
forrest@0	630	break;
forrest@0	631	case NC_CHAR:
forrest@0	632	if (!(val = (char )malloc(((len)+1) sizeof(char)))) {
forrest@0	633	perror("alloc_var: val");
forrest@0	634	exit(3);
forrest@0	635	}
forrest@0	636	break;
forrest@0	637	default:
forrest@0	638	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	639	exit(3);
forrest@0	640	}
forrest@0	641
forrest@0	642	return val;
forrest@0	643	}
forrest@0	644
forrest@0	645	void read_var(int ncid, int varid, nc_type nctype, int ndims, int dimids, struct dim **dim_idx)
forrest@0	646	{
forrest@0	647	int i;
forrest@0	648	size_t len = (size_t)1;
forrest@0	649	void *val;
forrest@0	650
forrest@0	651	/* Figure out the total size */
forrest@0	652	for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
forrest@0	653
forrest@0	654	val = alloc_var(nctype, len);
forrest@0	655
forrest@0	656	switch (nctype) {
forrest@0	657	case NC_FLOAT:
forrest@0	658	wrap_nc(nc_get_var_float(ncid, varid, val));
forrest@0	659	break;
forrest@0	660	case NC_INT:
forrest@0	661	wrap_nc(nc_get_var_int(ncid, varid, val));
forrest@0	662	break;
forrest@0	663	case NC_DOUBLE:
forrest@0	664	wrap_nc(nc_get_var_double(ncid, varid, val));
forrest@0	665	break;
forrest@0	666	case NC_CHAR:
forrest@0	667	wrap_nc(nc_get_var_text(ncid, varid, val));
forrest@0	668	break;
forrest@0	669	default:
forrest@0	670	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	671	exit(3);
forrest@0	672	}
forrest@0	673
forrest@0	674	return val;
forrest@0	675	}
forrest@0	676
forrest@0	677	void read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int dimids, struct dim **dim_idx, size_t tslice)
forrest@0	678	{
forrest@0	679	int i;
forrest@0	680	size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	681	void *val;
forrest@0	682
forrest@0	683	/* Make sure time is really the first dimension */
forrest@0	684	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
forrest@0	685	fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
forrest@0	686	exit(3);
forrest@0	687	}
forrest@0	688	/*
forrest@0	689	* Figure out the total size for the slice (start at second dimension)
forrest@0	690	* and build start/count vectors.
forrest@0	691	*/
forrest@0	692	start[0] = tslice;
forrest@0	693	count[0] = 1;
forrest@0	694	for (i = 1; i < ndims; i++) {
forrest@0	695	start[i] = 0;
forrest@0	696	count[i] = (dim_idx[dimids[i]])->len;
forrest@0	697	len *= (dim_idx[dimids[i]])->len;
forrest@0	698	}
forrest@0	699
forrest@0	700	val = alloc_var(nctype, len);
forrest@0	701
forrest@0	702	switch (nctype) {
forrest@0	703	case NC_FLOAT:
forrest@0	704	wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
forrest@0	705	break;
forrest@0	706	case NC_INT:
forrest@0	707	wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
forrest@0	708	break;
forrest@0	709	case NC_DOUBLE:
forrest@0	710	wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
forrest@0	711	break;
forrest@0	712	case NC_CHAR:
forrest@0	713	wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
forrest@0	714	break;
forrest@0	715	default:
forrest@0	716	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	717	exit(3);
forrest@0	718	}
forrest@0	719
forrest@0	720	return val;
forrest@0	721	}
forrest@0	722
forrest@0	723	void write_var(int ncid, int varid, nc_type nctype, void *val)
forrest@0	724	{
forrest@0	725	switch (nctype) {
forrest@0	726	case NC_FLOAT:
forrest@0	727	wrap_nc(nc_put_var_float(ncid, varid, val));
forrest@0	728	break;
forrest@0	729	case NC_INT:
forrest@0	730	wrap_nc(nc_put_var_int(ncid, varid, val));
forrest@0	731	break;
forrest@0	732	case NC_DOUBLE:
forrest@0	733	wrap_nc(nc_put_var_double(ncid, varid, val));
forrest@0	734	break;
forrest@0	735	case NC_CHAR:
forrest@0	736	wrap_nc(nc_put_var_text(ncid, varid, val));
forrest@0	737	break;
forrest@0	738	default:
forrest@0	739	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	740	exit(3);
forrest@0	741	}
forrest@0	742
forrest@0	743	return;
forrest@0	744	}
forrest@0	745
forrest@0	746	void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int dimids, struct dim dim_idx, void val, size_t tslice)
forrest@0	747	{
forrest@0	748	int i;
forrest@0	749	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	750
forrest@0	751	/* Make sure time is really the first dimension */
forrest@0	752	if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
forrest@0	753	fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
forrest@0	754	exit(3);
forrest@0	755	}
forrest@0	756
forrest@0	757	/* Build start/count vectors */
forrest@0	758	start[0] = tslice;
forrest@0	759	count[0] = 1;
forrest@0	760	for (i = 1; i < ndims; i++) {
forrest@0	761	start[i] = 0;
forrest@0	762	count[i] = (dim_idx[dimids[i]])->len;
forrest@0	763	}
forrest@0	764
forrest@0	765	switch (nctype) {
forrest@0	766	case NC_FLOAT:
forrest@0	767	wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
forrest@0	768	break;
forrest@0	769	case NC_INT:
forrest@0	770	wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
forrest@0	771	break;
forrest@0	772	case NC_DOUBLE:
forrest@0	773	wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
forrest@0	774	break;
forrest@0	775	case NC_CHAR:
forrest@0	776	wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
forrest@0	777	break;
forrest@0	778	default:
forrest@0	779	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	780	exit(3);
forrest@0	781	}
forrest@0	782
forrest@0	783	return;
forrest@0	784	}
forrest@0	785
forrest@0	786	void copy_metadata(int in_ncid, struct var *in_var_head,
forrest@1	787	struct dim *in_dim_idx, int out_ncid, struct var out_var_head,
forrest@1	788	double *tbounds)
forrest@0	789	{
forrest@0	790	int i, j;
forrest@0	791	float hr[HOURS_PER_DAY];
forrest@0	792	struct var in_vnode, out_vnode;
forrest@0	793	void *val;
forrest@0	794
forrest@0	795	for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
forrest@0	796
forrest@0	797	if (in_var_head) {
forrest@0	798	in_vnode = in_var_head;
forrest@0	799	/*
forrest@0	800	* March through input variables to stuff metadata values into
forrest@0	801	* the new output files. NOTE: Time-based metadata variables
forrest@0	802	* should contain only the last (time-slice) value (from all
forrest@0	803	* input files).
forrest@0	804	*/
forrest@0	805	for (j = 0; j == 0 \|\| in_vnode != in_var_head; j++) {
forrest@0	806	if (in_vnode->metadata) {
forrest@0	807	out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
forrest@0	808	if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
forrest@0	809	/* time is not the first dimension */
forrest@0	810	#ifdef DEBUG
forrest@0	811	printf("Copying metadata variable %s\n",
forrest@0	812	in_vnode->name);
forrest@0	813	#endif /* DEBUG */
forrest@0	814	val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
forrest@0	815	write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
forrest@0	816	}
forrest@0	817	else {
forrest@0	818	/* time is the first dimension */
forrest@0	819	#ifdef DEBUG
forrest@0	820	printf("Copying last value of \
forrest@0	821	time-based metadata variable %s\n", in_vnode->name);
forrest@0	822	#endif /* DEBUG */
forrest@0	823	val = read_timeslice(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, ((input_dim_idx[in_vnode->dimids[0]])->len - 1));
forrest@1	824	if (!strcmp(in_vnode->name, time_bounds_name))
forrest@1	825	write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, tbounds, 0);
forrest@1	826	else
forrest@1	827	write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
forrest@0	828	}
forrest@0	829	free(val);
forrest@0	830	/*
forrest@0	831	* Just after the "time" variable, write out
forrest@0	832	* the "hour" variable values.
forrest@0	833	*/
forrest@0	834	if (!strcmp(in_vnode->name, time_name)) {
forrest@0	835	out_vnode = varlist_find_by_name(out_var_head, hour_name);
forrest@0	836	write_var(out_ncid, out_vnode->ncvarid,
forrest@0	837	out_vnode->nctype, hr);
forrest@0	838	}
forrest@0	839	}
forrest@0	840	else {
forrest@0	841	#ifdef DEBUG
forrest@0	842	printf("Skipping variable %s\n",
forrest@0	843	in_vnode->name);
forrest@0	844	#endif /* DEBUG */
forrest@0	845	}
forrest@0	846	in_vnode = in_vnode->next;
forrest@0	847	}
forrest@0	848	}
forrest@0	849
forrest@0	850	return;
forrest@0	851	}
forrest@0	852
forrest@1	853	void get_time_bounds(char first_fname, char last_fname, double *tbounds)
forrest@1	854	{
forrest@1	855	int time_dimid, time_bounds_varid;
forrest@1	856	size_t time_len, time_bounds_index[2];
forrest@1	857	double bnd1, bnd2;
forrest@1	858
forrest@1	859	/* Open first input file */
forrest@1	860	wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
forrest@1	861	/* Get dimension ID of the time dimension */
forrest@1	862	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1	863	/* Get length of time dimension */
forrest@1	864	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1	865	/* Get variable ID of time_bounds variable */
forrest@1	866	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
forrest@1	867	/* Set index for reading the first value of time_bounds from the
forrest@1	868	* first timeslice of the first file */
forrest@1	869	time_bounds_index[0] = time_bounds_index[1] = 0;
forrest@1	870	/* Read the value */
forrest@1	871	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1	872	time_bounds_index, &bnd1));
forrest@1	873	/* If the first and last file are not the same, close the first one and
forrest@1	874	* open the second one */
forrest@1	875	if (strcmp(first_fname, last_fname)) {
forrest@1	876	/* Close the first input file */
forrest@1	877	wrap_nc(nc_close(input_ncid));
forrest@1	878	/* Open the last input file */
forrest@1	879	wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@1	880	/* Get dimension ID of the time dimension */
forrest@1	881	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1	882	/* Get length of time dimension */
forrest@1	883	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1	884	/* Get variable ID of time_bounds variable */
forrest@1	885	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
forrest@1	886	&time_bounds_varid));
forrest@1	887	}
forrest@1	888	/* Set index for reading the second value of time_bounds from the
forrest@1	889	* last timeslice of the last file */
forrest@1	890	time_bounds_index[0] = time_len - 1; time_bounds_index[1] = 1;
forrest@1	891	/* Read the value */
forrest@1	892	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1	893	time_bounds_index, &bnd2));
forrest@1	894
forrest@1	895	/* Close the last input file */
forrest@1	896	wrap_nc(nc_close(input_ncid));
forrest@1	897
forrest@1	898	tbounds[0] = bnd1;
forrest@1	899	tbounds[1] = bnd2;
forrest@1	900
forrest@1	901	return;
forrest@1	902	}
forrest@1	903
forrest@1	904	void open_inout(char first_fname, char last_fname, char mean_fname, char stddev_fname, char *flist, size_t nsamples)
forrest@0	905	{
forrest@0	906	char mean_history_gatt, stddev_history_gatt,
forrest@0	907	*mean_prefix="Statistical means from history files:",
forrest@0	908	*stddev_prefix="Statistical standard deviations from history files:";
forrest@1	909	double tbounds[2];
forrest@0	910
forrest@0	911	/*
forrest@0	912	* Construct strings for history global attributes for the two output
forrest@0	913	* files.
forrest@0	914	*/
forrest@0	915	if (!(mean_history_gatt = (char )malloc((strlen(mean_prefix) + strlen(flist)+1)sizeof(char)))) {
forrest@0	916	perror("open_inout:mean_history_gatt");
forrest@0	917	exit(2);
forrest@0	918	}
forrest@0	919	if (!(stddev_history_gatt = (char )malloc((strlen(stddev_prefix) + strlen(flist)+1)sizeof(char)))) {
forrest@0	920	perror("open_inout:stddev_history_gatt");
forrest@0	921	exit(2);
forrest@0	922	}
forrest@0	923	sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
forrest@0	924	sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
forrest@0	925
forrest@1	926	/* The two time_bounds values must be handled explicitly by obtaining
forrest@1	927	* the first one from the first file and the last one from the last
forrest@1	928	* file. These are then passed to copy_metadata() below. */
forrest@1	929	get_time_bounds(first_fname, last_fname, tbounds);
forrest@1	930	#ifdef DEBUG
forrest@1	931	fprintf(stderr, "Got back tbounds=%lf,%lf\n", tbounds[0], tbounds[1]);
forrest@1	932	#endif /* DEBUG */
forrest@1	933
forrest@1	934	/* Open last input file */
forrest@1	935	wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@0	936	/* Inquire about number of dimensions, variables, global attributes
forrest@0	937	* and the ID of the unlimited dimension
forrest@0	938	*/
forrest@0	939	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
forrest@0	940	&input_unlimdimid));
forrest@0	941	/* Grab dimension IDs and lengths from input file */
forrest@0	942	get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
forrest@0	943	/* Grab desired global attributes from input file */
forrest@0	944	get_gatts(input_ncid, input_ngatts, &input_gatt_head);
forrest@0	945	/* Grab variable IDs and attributes from input file */
forrest@0	946	get_vars(input_ncid, input_nvars, &input_var_head);
forrest@0	947	#ifdef DEBUG
forrest@0	948	varlist_print(input_var_head);
forrest@0	949	#endif /* DEBUG */
forrest@0	950	/* Create netCDF files for monthly means and standard deviations */
forrest@0	951	/* Create new netCDF files */
forrest@0	952	wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
forrest@0	953	wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
forrest@0	954	/* Define dimensions */
forrest@0	955	mean_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0	956	input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
forrest@0	957	&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
forrest@0	958	stddev_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0	959	input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
forrest@0	960	&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
forrest@0	961	/* Define variables and their attributes */
forrest@0	962	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0	963	idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
forrest@0	964	MEAN_TYPE);
forrest@0	965	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0	966	idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
forrest@0	967	STDDEV_TYPE);
forrest@0	968	/* Store global attributes */
forrest@0	969	mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
forrest@0	970	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
forrest@0	971	stddev_history_gatt);
forrest@0	972	/* End define mode */
forrest@0	973	wrap_nc(nc_enddef(mean_ncid));
forrest@0	974	wrap_nc(nc_enddef(stddev_ncid));
forrest@0	975	/* Write out metdata variables that do not depend on "time" */
forrest@0	976	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
forrest@1	977	mean_var_head, tbounds);
forrest@0	978	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
forrest@1	979	stddev_var_head, tbounds);
forrest@0	980
forrest@0	981	wrap_nc(nc_close(input_ncid));
forrest@0	982
forrest@0	983	return;
forrest@0	984	}
forrest@0	985
forrest@0	986	size_t count_samples(int nifnames, char **input_fnames)
forrest@0	987	{
forrest@0	988	int i;
forrest@0	989	char name[NC_MAX_NAME+1];
forrest@0	990	size_t len, cnt = 0;
forrest@0	991
forrest@0	992	for (i = 0; i < nifnames; i++) {
forrest@0	993	/* Open input file */
forrest@0	994	wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0	995	/*
forrest@0	996	* Inquire about number of dimensions, variables, global
forrest@0	997	* attributes and the ID of the unlimited dimension
forrest@0	998	*/
forrest@0	999	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0	1000	&input_ngatts, &input_unlimdimid));
forrest@0	1001	wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0	1002	if (strcmp(name, time_name)) {
forrest@0	1003	fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0	1004	exit(4);
forrest@0	1005	}
forrest@0	1006	#ifdef DEBUG
forrest@0	1007	printf("%ld samples in %s\n", (long)len, input_fnames[i]);
forrest@0	1008	#endif /* DEBUG */
forrest@0	1009	wrap_nc(nc_close(input_ncid));
forrest@0	1010	cnt += len;
forrest@0	1011	}
forrest@0	1012	return cnt;
forrest@0	1013	}
forrest@0	1014
forrest@0	1015	void alloc_means_stddevs(size_t d1, size_t d2, double *meansp, double stddevsp, size_t **cell_samplesp)
forrest@0	1016	{
forrest@0	1017	/*
forrest@0	1018	* Allocate space for arrays of means and standard deviations by
forrest@0	1019	* hour of day.
forrest@0	1020	*/
forrest@0	1021	int i;
forrest@0	1022	size_t **cell_samples;
forrest@0	1023	double means, stddevs;
forrest@0	1024
forrest@0	1025	if (!(meansp = (double )malloc(HOURS_PER_DAY sizeof(double *)))) {
forrest@0	1026	perror("alloc_means_stddevs:*meansp");
forrest@0	1027	exit(2);
forrest@0	1028	}
forrest@0	1029	if (!(stddevsp = (double )malloc(HOURS_PER_DAY sizeof(double *)))) {
forrest@0	1030	perror("alloc_means_stddevs:*stddevsp");
forrest@0	1031	exit(2);
forrest@0	1032	}
forrest@0	1033	if (!(cell_samplesp = (size_t )malloc(HOURS_PER_DAY sizeof(size_t *)))) {
forrest@0	1034	perror("alloc_means_stddevs:*cell_samplesp");
forrest@0	1035	exit(2);
forrest@0	1036	}
forrest@0	1037
forrest@0	1038	means = *meansp;
forrest@0	1039	stddevs = *stddevsp;
forrest@0	1040	cell_samples = *cell_samplesp;
forrest@0	1041
forrest@0	1042	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0	1043	if (!(means[i] = (double )malloc(d1 d2 * sizeof(double)))) {
forrest@0	1044	perror("alloc_means_stddevs:means[i]");
forrest@0	1045	exit(2);
forrest@0	1046	}
forrest@0	1047	if (!(stddevs[i] = (double )malloc(d1 d2 * sizeof(double)))) {
forrest@0	1048	perror("alloc_means_stddevs:stddevs[i]");
forrest@0	1049	exit(2);
forrest@0	1050	}
forrest@0	1051	if (!(cell_samples[i] = (size_t )malloc(d1 d2 * sizeof(size_t)))) {
forrest@0	1052	perror("alloc_means_stddevs:cell_samples[i]");
forrest@0	1053	exit(2);
forrest@0	1054	}
forrest@0	1055	}
forrest@0	1056
forrest@0	1057	return;
forrest@0	1058	}
forrest@0	1059
forrest@0	1060	void init_means_stddevs(size_t d1, size_t d2, double means, double stddevs, size_t **cell_samples, float FillValue)
forrest@0	1061	{
forrest@0	1062	int hr, i, j, pos;
forrest@0	1063
forrest@0	1064	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1065	for (i = 0; i < d1; i++) {
forrest@0	1066	#pragma _CRI concurrent
forrest@0	1067	for (j = 0; j < d2; j++) {
forrest@0	1068	pos = i*d2+j;
forrest@0	1069	means[hr][pos] = FillValue;
forrest@0	1070	stddevs[hr][pos] = FillValue;
forrest@0	1071	cell_samples[hr][pos] = 0;
forrest@0	1072	}
forrest@0	1073	}
forrest@0	1074	}
forrest@0	1075
forrest@0	1076	return;
forrest@0	1077	}
forrest@0	1078
forrest@0	1079	void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
forrest@0	1080	{
forrest@0	1081	int i, j, hr, pos;
forrest@0	1082
forrest@0	1083	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1084	for (i = 0; i < d1; i++) {
forrest@0	1085	#pragma _CRI concurrent
forrest@0	1086	for (j = 0; j < d2; j++) {
forrest@0	1087	pos = i*d2+j;
forrest@0	1088	cell_samples[hr][pos] = 0;
forrest@0	1089	}
forrest@0	1090	}
forrest@0	1091	}
forrest@0	1092
forrest@0	1093	return;
forrest@0	1094	}
forrest@0	1095
forrest@0	1096	void free_means_stddevs(double means, double stddevs, size_t **cell_samples)
forrest@0	1097	{
forrest@0	1098	/*
forrest@0	1099	* Free space from arrays of means and standard deviations by
forrest@0	1100	* hour of day.
forrest@0	1101	*/
forrest@0	1102	int i;
forrest@0	1103
forrest@0	1104	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0	1105	free(means[i]);
forrest@0	1106	free(stddevs[i]);
forrest@0	1107	free(cell_samples[i]);
forrest@0	1108	}
forrest@0	1109
forrest@0	1110	free(means);
forrest@0	1111	free(stddevs);
forrest@0	1112	free(cell_samples);
forrest@0	1113
forrest@0	1114	return;
forrest@0	1115	}
forrest@0	1116
forrest@0	1117	void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
forrest@0	1118	char FillFlag, float FillValue, double means, size_t cell_samples)
forrest@0	1119	{
forrest@0	1120	int i, j, hr, pos;
forrest@0	1121
forrest@0	1122	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1123
forrest@0	1124	for (i = 0; i < d1; i++) {
forrest@0	1125	#pragma _CRI concurrent
forrest@0	1126	for (j = 0; j < d2; j++) {
forrest@0	1127	pos = i*d2+j;
forrest@0	1128	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue)) {
forrest@0	1129	if (cell_samples[hr][pos] == 0)
forrest@0	1130	means[hr][pos] = (double)val[pos];
forrest@0	1131	else
forrest@0	1132	means[hr][pos] += (double)val[pos];
forrest@0	1133	++cell_samples[hr][pos];
forrest@0	1134	}
forrest@0	1135	}
forrest@0	1136	}
forrest@0	1137
forrest@0	1138	return;
forrest@0	1139	}
forrest@0	1140
forrest@0	1141	void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
forrest@0	1142	char FillFlag, float FillValue, double means, size_t cell_samples)
forrest@0	1143	{
forrest@0	1144	int i, j, hr, pos;
forrest@0	1145
forrest@0	1146	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1147
forrest@0	1148	for (i = 0; i < d1; i++) {
forrest@0	1149	#pragma _CRI concurrent
forrest@0	1150	for (j = 0; j < d2; j++) {
forrest@0	1151	pos = i*d2+j;
forrest@0	1152	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue)) {
forrest@0	1153	if (cell_samples[hr][pos] == 0)
forrest@0	1154	means[hr][pos] = val[pos];
forrest@0	1155	else
forrest@0	1156	means[hr][pos] += val[pos];
forrest@0	1157	++cell_samples[hr][pos];
forrest@0	1158	}
forrest@0	1159	}
forrest@0	1160	}
forrest@0	1161
forrest@0	1162	return;
forrest@0	1163	}
forrest@0	1164
forrest@0	1165
forrest@0	1166	void divide_means(size_t d1, size_t d2, double means, size_t cell_samples)
forrest@0	1167	{
forrest@0	1168	int i, j, hr, pos;
forrest@0	1169
forrest@0	1170	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1171	for (i = 0; i < d1; i++) {
forrest@0	1172	#pragma _CRI concurrent
forrest@0	1173	for (j = 0; j < d2; j++) {
forrest@0	1174	pos = i*d2+j;
forrest@0	1175	if (cell_samples[hr][pos] != 0) {
forrest@0	1176	means[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0	1177	}
forrest@0	1178	}
forrest@0	1179	}
forrest@0	1180	}
forrest@0	1181
forrest@0	1182	return;
forrest@0	1183	}
forrest@0	1184
forrest@0	1185	void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
forrest@0	1186	char FillFlag, float FillValue, double means, double stddevs,
forrest@0	1187	size_t **cell_samples)
forrest@0	1188	{
forrest@0	1189	int i, j, hr, pos;
forrest@0	1190	double delta;
forrest@0	1191
forrest@0	1192	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1193
forrest@0	1194	for (i = 0; i < d1; i++) {
forrest@0	1195	#pragma _CRI concurrent
forrest@0	1196	for (j = 0; j < d2; j++) {
forrest@0	1197	pos = i*d2+j;
forrest@0	1198	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue
forrest@0	1199	&& means[hr][pos] != FillValue)) {
forrest@0	1200	delta = means[hr][pos] - (double)val[pos];
forrest@0	1201	if (cell_samples[hr][pos] == 0)
forrest@0	1202	stddevs[hr][pos] = delta * delta;
forrest@0	1203	else
forrest@0	1204	stddevs[hr][pos] += delta * delta;
forrest@0	1205	++cell_samples[hr][pos];
forrest@0	1206	}
forrest@0	1207	}
forrest@0	1208	}
forrest@0	1209
forrest@0	1210	return;
forrest@0	1211	}
forrest@0	1212
forrest@0	1213	void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
forrest@0	1214	char FillFlag, float FillValue, double means, double stddevs,
forrest@0	1215	size_t **cell_samples)
forrest@0	1216	{
forrest@0	1217	int i, j, hr, pos;
forrest@0	1218	double delta;
forrest@0	1219
forrest@0	1220	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1221
forrest@0	1222	for (i = 0; i < d1; i++) {
forrest@0	1223	#pragma _CRI concurrent
forrest@0	1224	for (j = 0; j < d2; j++) {
forrest@0	1225	pos = i*d2+j;
forrest@0	1226	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue
forrest@0	1227	&& means[hr][pos] != FillValue)) {
forrest@0	1228	delta = means[hr][pos] - val[pos];
forrest@0	1229	if (cell_samples[hr][pos] == 0)
forrest@0	1230	stddevs[hr][pos] = delta * delta;
forrest@0	1231	else
forrest@0	1232	stddevs[hr][pos] += delta * delta;
forrest@0	1233	++cell_samples[hr][pos];
forrest@0	1234	}
forrest@0	1235	}
forrest@0	1236	}
forrest@0	1237
forrest@0	1238	return;
forrest@0	1239	}
forrest@0	1240
forrest@0	1241
forrest@0	1242	void divide_sqrt_stddevs(size_t d1, size_t d2, double stddevs, size_t cell_samples)
forrest@0	1243	{
forrest@0	1244	int i, j, hr, pos;
forrest@0	1245
forrest@0	1246	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1247	for (i = 0; i < d1; i++) {
forrest@0	1248	#pragma _CRI concurrent
forrest@0	1249	for (j = 0; j < d2; j++) {
forrest@0	1250	pos = i*d2+j;
forrest@0	1251	if (cell_samples[hr][pos] != 0) {
forrest@0	1252	stddevs[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0	1253	stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
forrest@0	1254	}
forrest@0	1255	}
forrest@0	1256	}
forrest@0	1257	}
forrest@0	1258
forrest@0	1259	return;
forrest@0	1260	}
forrest@0	1261
forrest@0	1262	void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void *valsp, int mcsecp)
forrest@0	1263	{
forrest@0	1264	int i;
forrest@0	1265	void **vals;
forrest@0	1266
forrest@0	1267	/* Allocate space for all timeslices */
forrest@0	1268	if (!(valsp = (void )malloc(nsamples sizeof(void *)))) {
forrest@0	1269	perror("alloc_all_samples:*valsp");
forrest@0	1270	exit(2);
forrest@0	1271	}
forrest@0	1272	vals = *valsp;
forrest@0	1273	for (i = 0; i < nsamples; i++) {
forrest@0	1274	vals[i] = alloc_var(nctype, (d1 * d2));
forrest@0	1275	}
forrest@0	1276	if (!(mcsecp = (int )malloc(nsamples * sizeof(int)))) {
forrest@0	1277	perror("alloc_all_samples:*mcsecp");
forrest@0	1278	exit(2);
forrest@0	1279	}
forrest@0	1280
forrest@0	1281	return;
forrest@0	1282	}
forrest@0	1283
forrest@0	1284	void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void *vals, int mcsec)
forrest@0	1285	{
forrest@0	1286	int i, j;
forrest@0	1287	float **fvals = NULL;
forrest@0	1288	double **dvals = NULL;
forrest@0	1289
forrest@0	1290	switch(nctype) {
forrest@0	1291	case NC_FLOAT:
forrest@0	1292	fvals = (float **)vals;
forrest@0	1293	break;
forrest@0	1294	case NC_DOUBLE:
forrest@0	1295	dvals = (double **)vals;
forrest@0	1296	break;
forrest@0	1297	default:
forrest@0	1298	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1299	exit(3);
forrest@0	1300	}
forrest@0	1301
forrest@0	1302	for (i = 0; i < nsamples; i++) {
forrest@0	1303	for (j = 0; j < (d1 * d2); j++) {
forrest@0	1304	switch(nctype) {
forrest@0	1305	case NC_FLOAT:
forrest@0	1306	fvals[i][j] = (float)0.0;
forrest@0	1307	break;
forrest@0	1308	case NC_DOUBLE:
forrest@0	1309	dvals[i][j] = (double)0.0;
forrest@0	1310	break;
forrest@0	1311	default:
forrest@0	1312	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1313	exit(3);
forrest@0	1314	}
forrest@0	1315	}
forrest@0	1316	mcsec[i] = 0;
forrest@0	1317	}
forrest@0	1318
forrest@0	1319	return;
forrest@0	1320	}
forrest@0	1321
forrest@0	1322	void free_all_samples(size_t nsamples, void *vals, int mcsec)
forrest@0	1323	{
forrest@0	1324	int i;
forrest@0	1325
forrest@0	1326	for (i = 0; i < nsamples; i++)
forrest@0	1327	if (vals[i]) free(vals[i]);
forrest@0	1328
forrest@0	1329	free(vals);
forrest@0	1330	free(mcsec);
forrest@0	1331
forrest@0	1332	return;
forrest@0	1333	}
forrest@0	1334
forrest@0	1335	void read_all_samples(int nifnames, char *input_fnames, size_t d1, size_t d2, size_t nsamples, char var_name, nc_type nctype, int level, int ndims, void *vals, int mcsec)
forrest@0	1336	{
forrest@0	1337	int i, ts;
forrest@0	1338	int varid;
forrest@0	1339	size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	1340	char name[NC_MAX_NAME+1];
forrest@0	1341
forrest@0	1342	for (i = 0; i < nifnames; i++) {
forrest@0	1343	#ifdef DEBUG
forrest@0	1344	printf("\tOpening %s", input_fnames[i]);
forrest@0	1345	if (ndims > 3) printf(" to retrieve level %d\n", level);
forrest@0	1346	#endif /* DEBUG */
forrest@0	1347	/* Open input file */
forrest@0	1348	wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0	1349	/*
forrest@0	1350	* Inquire about number of dimensions, variables, global
forrest@0	1351	* attributes and the ID of the unlimited dimension
forrest@0	1352	*/
forrest@0	1353	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0	1354	&input_ngatts, &input_unlimdimid));
forrest@0	1355	wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0	1356	if (strcmp(name, time_name)) {
forrest@0	1357	fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0	1358	exit(4);
forrest@0	1359	}
forrest@0	1360	/* Make sure we don't run off the end of allocated space */
forrest@0	1361	if ((nslice+len) > nsamples) {
forrest@0	1362	fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);
forrest@0	1363	exit(4);
forrest@0	1364	}
forrest@0	1365	/* Get seconds of day variable */
forrest@0	1366	wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
forrest@0	1367	wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));
forrest@0	1368	/* Get variable ID for requested variable */
forrest@0	1369	wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
forrest@0	1370	/* Retrieve time slice of desired variable */
forrest@0	1371	for (ts = 0; ts < len; ts++) {
forrest@0	1372	start[0] = ts;
forrest@0	1373	count[0] = 1;
forrest@0	1374	start[(ndims-2)] = 0;
forrest@0	1375	count[(ndims-2)] = d1;
forrest@0	1376	start[(ndims-1)] = 0;
forrest@0	1377	count[(ndims-1)] = d2;
forrest@0	1378	if (ndims == 4) {
forrest@0	1379	start[1] = level;
forrest@0	1380	count[1] = 1;
forrest@0	1381	}
forrest@0	1382	switch(nctype) {
forrest@0	1383	case NC_FLOAT:
forrest@0	1384	wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));
forrest@0	1385	break;
forrest@0	1386	case NC_DOUBLE:
forrest@0	1387	wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));
forrest@0	1388	break;
forrest@0	1389	default:
forrest@0	1390	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1391	exit(3);
forrest@0	1392	}
forrest@0	1393	nslice++;
forrest@0	1394	}
forrest@0	1395
forrest@0	1396	/* Close input file */
forrest@0	1397	wrap_nc(nc_close(input_ncid));
forrest@0	1398	}
forrest@0	1399
forrest@0	1400	return;
forrest@0	1401	}
forrest@0	1402
forrest@0	1403	void compute(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void *vals, int mcsec, char FillFlag, float FillValue, int stat_type, double means, double stddevs, size_t **cell_samples)
forrest@0	1404	{
forrest@0	1405	int i;
forrest@0	1406
forrest@0	1407	for (i = 0; i < nsamples; i++) {
forrest@0	1408	switch(nctype) {
forrest@0	1409	case NC_FLOAT:
forrest@0	1410	if (stat_type == MEAN_TYPE)
forrest@0	1411	add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0	1412	else
forrest@0	1413	add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0	1414	break;
forrest@0	1415	case NC_DOUBLE:
forrest@0	1416	if (stat_type == MEAN_TYPE)
forrest@0	1417	add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0	1418	else
forrest@0	1419	add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0	1420	break;
forrest@0	1421	default:
forrest@0	1422	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1423	exit(3);
forrest@0	1424	}
forrest@0	1425	}
forrest@0	1426
forrest@0	1427	/* Divide sums by number of samples */
forrest@0	1428	if (stat_type == MEAN_TYPE)
forrest@0	1429	divide_means(d1, d2, means, cell_samples);
forrest@0	1430	else
forrest@0	1431	divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
forrest@0	1432
forrest@0	1433	return;
forrest@0	1434	}
forrest@0	1435
forrest@0	1436	float double_to_float(size_t len, double dvar)
forrest@0	1437	{
forrest@0	1438	int i;
forrest@0	1439	float *fvar;
forrest@0	1440
forrest@0	1441	if (!(fvar = (float )malloc(len sizeof(float)))) {
forrest@0	1442	perror("double_to_float:fvar");
forrest@0	1443	exit(2);
forrest@0	1444	}
forrest@0	1445
forrest@0	1446	for (i = 0; i < len; i++)
forrest@0	1447	fvar[i] = (float)dvar[i];
forrest@0	1448
forrest@0	1449	return fvar;
forrest@0	1450	}
forrest@0	1451
forrest@0	1452	void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
forrest@0	1453	int level, int ndims, double **var_hours)
forrest@0	1454	{
forrest@0	1455	/* Output dimensions should be one larger than input dimensions */
forrest@0	1456	int i, hr;
forrest@0	1457	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	1458	float *fvar = NULL;
forrest@0	1459
forrest@0	1460	if (nctype == NC_FLOAT) {
forrest@0	1461	if (!(fvar = (float )malloc(d1 d2 * sizeof(float)))) {
forrest@0	1462	perror("write_var_hours:fvar");
forrest@0	1463	exit(2);
forrest@0	1464	}
forrest@0	1465	}
forrest@0	1466
forrest@0	1467	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1468	start[0] = 0;
forrest@0	1469	count[0] = 1; /* time */
forrest@0	1470	start[1] = hr;
forrest@0	1471	count[1] = 1; /* hour */
forrest@0	1472	start[(ndims-2)] = 0;
forrest@0	1473	count[(ndims-2)] = d1;
forrest@0	1474	start[(ndims-1)] = 0;
forrest@0	1475	count[(ndims-1)] = d2;
forrest@0	1476	if (ndims == 5) {
forrest@0	1477	start[2] = level;
forrest@0	1478	count[2] = 1;
forrest@0	1479	}
forrest@0	1480	switch (nctype) {
forrest@0	1481	case NC_FLOAT:
forrest@0	1482	for (i = 0; i < (d1 * d2); i++)
forrest@0	1483	fvar[i] = (float)var_hours[hr][i];
forrest@0	1484	wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
forrest@0	1485	break;
forrest@0	1486	case NC_DOUBLE:
forrest@0	1487	wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
forrest@0	1488	break;
forrest@0	1489	default:
forrest@0	1490	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1491	exit(3);
forrest@0	1492	}
forrest@0	1493	}
forrest@0	1494
forrest@0	1495	if (nctype == NC_FLOAT)
forrest@0	1496	free(fvar);
forrest@0	1497
forrest@0	1498	return;
forrest@0	1499	}
forrest@0	1500
forrest@0	1501	void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
forrest@0	1502	{
forrest@0	1503	int j, k, nlevels, *mcsec;
forrest@0	1504	size_t d1, d2, **cell_samples;
forrest@0	1505	double means, stddevs;
forrest@0	1506	struct var in_vnode, out_vnode;
forrest@0	1507	void **vals;
forrest@0	1508
forrest@0	1509	if (input_var_head) {
forrest@0	1510	in_vnode = input_var_head;
forrest@0	1511	/* March through non-metadata variables to compute stats */
forrest@0	1512	for (j = 0; j == 0 \|\| in_vnode != input_var_head; j++) {
forrest@0	1513	if (!in_vnode->metadata) {
forrest@0	1514	/* Make sure time is really the first dimension */
forrest@0	1515	if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
forrest@0	1516	fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
forrest@0	1517	exit(3);
forrest@0	1518	}
forrest@0	1519	/* Figure out input dimensions */
forrest@0	1520	if (in_vnode->ndims < 3 \|\| in_vnode->ndims > 4) {
forrest@0	1521	fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
forrest@0	1522	exit(3);
forrest@0	1523	}
forrest@0	1524	/* Assume 2D output; find dimensions */
forrest@0	1525	d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
forrest@0	1526	d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
forrest@0	1527	if (in_vnode->ndims == 3)
forrest@0	1528	nlevels = 1;
forrest@0	1529	else
forrest@0	1530	nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
forrest@0	1531	/* Allocate working space for means and stddevs */
forrest@0	1532	alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
forrest@0	1533	init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
forrest@0	1534	/* Allocate and initize space for entire field across all files */
forrest@0	1535	alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);
forrest@0	1536	init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);
forrest@0	1537	printf("Computing statistics for %s\n",
forrest@0	1538	in_vnode->name);
forrest@0	1539	/* Compute means and stddevs, then write them */
forrest@0	1540	for (k = 0; k < nlevels; k++) {
forrest@0	1541	/* Read the entire fields from all files */
forrest@0	1542	read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);
forrest@0	1543	/* Compute means */
forrest@0	1544	compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
forrest@0	1545	/* Find corresponding output variable on the mean output file */
forrest@0	1546	out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
forrest@0	1547	/* Write out the means for this variable */
forrest@0	1548	write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
forrest@0	1549	/* Zero out cell_samples so they can be used as a flag again for computing stddevs */
forrest@0	1550	reset_cell_samples(d1, d2, cell_samples);
forrest@0	1551	/* Compute stddevs using means */
forrest@0	1552	compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
forrest@0	1553	/* Find corresponding output variable on the stddev output file */
forrest@0	1554	out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
forrest@0	1555	/* Write out stddevs for this variable */
forrest@0	1556	write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
forrest@0	1557	}
forrest@0	1558
forrest@0	1559	/* Free all samples */
forrest@0	1560	free_all_samples(nsamples, vals, mcsec);
forrest@0	1561	/* Free working space for means and stddevs */
forrest@0	1562	free_means_stddevs(means, stddevs, cell_samples);
forrest@0	1563	}
forrest@0	1564	in_vnode = in_vnode->next;
forrest@0	1565	}
forrest@0	1566	}
forrest@0	1567	return;
forrest@0	1568	}
forrest@0	1569
forrest@0	1570	void usage(char *program)
forrest@0	1571	{
forrest@0	1572	fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
forrest@0	1573	fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
forrest@0	1574	exit(4);
forrest@0	1575	}
forrest@0	1576
forrest@0	1577	int main(int argc, char **argv)
forrest@0	1578	{
forrest@0	1579	int i, ic, nifnames;
forrest@0	1580	size_t len, pos, nsamples;
forrest@0	1581	char mfname, sfname, *ifnames, flist;
forrest@0	1582
forrest@0	1583	ifnames = NULL;
forrest@0	1584	mfname = sfname = flist = NULL;
forrest@0	1585	nifnames = 0;
forrest@0	1586
forrest@0	1587	#ifdef DEBUG
forrest@0	1588	printf("Number of metadata variables (nmvars) = %d\n", nmvars);
forrest@0	1589	fflush(stdout);
forrest@0	1590	#endif /* DEBUG */
forrest@0	1591
forrest@0	1592	/* check command-line flags and switches */
forrest@0	1593	while ((ic = getopt(argc, argv, "m:s:")) != -1) {
forrest@0	1594	switch(ic) {
forrest@0	1595	case 'm':
forrest@0	1596	if (!(mfname = strdup(optarg))) {
forrest@0	1597	perror("mfname");
forrest@0	1598	exit(2);
forrest@0	1599	}
forrest@0	1600	break;
forrest@0	1601	case 's':
forrest@0	1602	if (!(sfname = strdup(optarg))) {
forrest@0	1603	perror("sfname");
forrest@0	1604	exit(2);
forrest@0	1605	}
forrest@0	1606	break;
forrest@0	1607	}
forrest@0	1608	}
forrest@0	1609
forrest@0	1610	if (!mfname) {
forrest@0	1611	fprintf(stderr, "Output file name for writing means required.\n");
forrest@0	1612	usage(argv[0]);
forrest@0	1613	}
forrest@0	1614	if (!sfname) {
forrest@0	1615	fprintf(stderr, "Output file name for writing standard deviations required.\n");
forrest@0	1616	usage(argv[0]);
forrest@0	1617	}
forrest@0	1618
forrest@0	1619	if (optind < argc) {
forrest@0	1620	if (!(ifnames = (char *)malloc((argc-optind+1)sizeof(char *)))) {
forrest@0	1621	perror("ifnames");
forrest@0	1622	exit(2);
forrest@0	1623	}
forrest@0	1624	for (i = optind; i < argc; i++) {
forrest@0	1625	if (!(ifnames[nifnames++] = strdup(argv[i]))) {
forrest@0	1626	perror("ifnames[nifnames++]");
forrest@0	1627	exit(2);
forrest@0	1628	}
forrest@0	1629	}
forrest@0	1630	}
forrest@0	1631	else {
forrest@0	1632	fprintf(stderr, "At least one input file name is required.\n");
forrest@0	1633	usage(argv[0]);
forrest@0	1634	}
forrest@0	1635
forrest@0	1636
forrest@0	1637	/*
forrest@0	1638	* Build list of input files to be included in the global history
forrest@0	1639	* attribute on the two outputs files.
forrest@0	1640	*/
forrest@0	1641	for (i = len = 0; i < nifnames; i++)
forrest@0	1642	len += strlen(ifnames[i]);
forrest@0	1643	len += nifnames + 1; /* space in front of every name + null terminator */
forrest@0	1644	if (!(flist = (char )malloc(len sizeof(char)))) {
forrest@0	1645	perror("flist");
forrest@0	1646	exit(2);
forrest@0	1647	}
forrest@0	1648	for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
forrest@0	1649	sprintf(&flist[pos], " %s", ifnames[i]);
forrest@0	1650	#ifdef DEBUG
forrest@0	1651	printf("flist=%s\n", flist);
forrest@0	1652	fflush(stdout);
forrest@0	1653	#endif /* DEBUG */
forrest@0	1654
forrest@0	1655	/* Open every file to count up the number of time samples. */
forrest@0	1656	nsamples = count_samples(nifnames, ifnames);
forrest@0	1657	#ifdef DEBUG
forrest@0	1658	printf("Number of samples across all files = %ld\n", (long)nsamples);
forrest@0	1659	fflush(stdout);
forrest@0	1660	#endif /* DEBUG */
forrest@0	1661
forrest@0	1662	/*
forrest@0	1663	* Open the last input file and the two output files (for mean and
forrest@0	1664	* standard deviation). Define dimensions, variables, and attributes
forrest@0	1665	* for the two output files based on the last input files. The
forrest@0	1666	* last file is used because some metadata variables will contain
forrest@0	1667	* only values from the last time slice from the period over which
forrest@0	1668	* the statistics are computed.
forrest@0	1669	*/
forrest@1	1670	open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
forrest@1	1671	nsamples);
forrest@0	1672
forrest@0	1673	compute_stats(nifnames, ifnames, nsamples);
forrest@0	1674
forrest@0	1675	/* Close the two output files */
forrest@0	1676	wrap_nc(nc_close(mean_ncid));
forrest@0	1677	wrap_nc(nc_close(stddev_ncid));
forrest@0	1678
forrest@0	1679	return 0;
forrest@0	1680	}

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 01 Oct 2007 15:12:14 -0400
changeset 1	2ce4ee911439
parent 0	3c02cce30be8
child 3	d3122367777b
permissions	-rw-r--r--