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

author	Forrest Hoffman <forrest@climatemodeling.org>
	Wed, 26 Sep 2007 17:16:40 -0400
changeset 0	3c02cce30be8
child 1	2ce4ee911439
permissions	-rw-r--r--