c-lamp_postproc: h1_summary2.c@dd8e6719647b (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@3	791	float hr[HOURS_PER_DAY], mean_time;
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@3	826	else if (!strcmp(in_vnode->name, time_name)) {
forrest@3	827	/* force the timestamp to be the
forrest@3	828	* mean of the time bounds */
forrest@3	829	mean_time = (tbounds[0] + tbounds[1]) / 2.0;
forrest@3	830	write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, &mean_time, 0);
forrest@3	831	}
forrest@1	832	else
forrest@1	833	write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
forrest@0	834	}
forrest@0	835	free(val);
forrest@0	836	/*
forrest@0	837	* Just after the "time" variable, write out
forrest@0	838	* the "hour" variable values.
forrest@0	839	*/
forrest@0	840	if (!strcmp(in_vnode->name, time_name)) {
forrest@0	841	out_vnode = varlist_find_by_name(out_var_head, hour_name);
forrest@0	842	write_var(out_ncid, out_vnode->ncvarid,
forrest@0	843	out_vnode->nctype, hr);
forrest@0	844	}
forrest@0	845	}
forrest@0	846	else {
forrest@0	847	#ifdef DEBUG
forrest@0	848	printf("Skipping variable %s\n",
forrest@0	849	in_vnode->name);
forrest@0	850	#endif /* DEBUG */
forrest@0	851	}
forrest@0	852	in_vnode = in_vnode->next;
forrest@0	853	}
forrest@0	854	}
forrest@0	855
forrest@0	856	return;
forrest@0	857	}
forrest@0	858
forrest@1	859	void get_time_bounds(char first_fname, char last_fname, double *tbounds)
forrest@1	860	{
forrest@1	861	int time_dimid, time_bounds_varid;
forrest@1	862	size_t time_len, time_bounds_index[2];
forrest@1	863	double bnd1, bnd2;
forrest@1	864
forrest@1	865	/* Open first input file */
forrest@1	866	wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
forrest@1	867	/* Get dimension ID of the time dimension */
forrest@1	868	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1	869	/* Get length of time dimension */
forrest@1	870	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1	871	/* Get variable ID of time_bounds variable */
forrest@1	872	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
forrest@1	873	/* Set index for reading the first value of time_bounds from the
forrest@1	874	* first timeslice of the first file */
forrest@1	875	time_bounds_index[0] = time_bounds_index[1] = 0;
forrest@1	876	/* Read the value */
forrest@1	877	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1	878	time_bounds_index, &bnd1));
forrest@1	879	/* If the first and last file are not the same, close the first one and
forrest@1	880	* open the second one */
forrest@1	881	if (strcmp(first_fname, last_fname)) {
forrest@1	882	/* Close the first input file */
forrest@1	883	wrap_nc(nc_close(input_ncid));
forrest@1	884	/* Open the last input file */
forrest@1	885	wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@1	886	/* Get dimension ID of the time dimension */
forrest@1	887	wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
forrest@1	888	/* Get length of time dimension */
forrest@1	889	wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
forrest@1	890	/* Get variable ID of time_bounds variable */
forrest@1	891	wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
forrest@1	892	&time_bounds_varid));
forrest@1	893	}
forrest@1	894	/* Set index for reading the second value of time_bounds from the
forrest@1	895	* last timeslice of the last file */
forrest@1	896	time_bounds_index[0] = time_len - 1; time_bounds_index[1] = 1;
forrest@1	897	/* Read the value */
forrest@1	898	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,
forrest@1	899	time_bounds_index, &bnd2));
forrest@1	900
forrest@1	901	/* Close the last input file */
forrest@1	902	wrap_nc(nc_close(input_ncid));
forrest@1	903
forrest@1	904	tbounds[0] = bnd1;
forrest@1	905	tbounds[1] = bnd2;
forrest@1	906
forrest@1	907	return;
forrest@1	908	}
forrest@1	909
forrest@1	910	void open_inout(char first_fname, char last_fname, char mean_fname, char stddev_fname, char *flist, size_t nsamples)
forrest@0	911	{
forrest@0	912	char mean_history_gatt, stddev_history_gatt,
forrest@0	913	*mean_prefix="Statistical means from history files:",
forrest@0	914	*stddev_prefix="Statistical standard deviations from history files:";
forrest@1	915	double tbounds[2];
forrest@0	916
forrest@0	917	/*
forrest@0	918	* Construct strings for history global attributes for the two output
forrest@0	919	* files.
forrest@0	920	*/
forrest@0	921	if (!(mean_history_gatt = (char )malloc((strlen(mean_prefix) + strlen(flist)+1)sizeof(char)))) {
forrest@0	922	perror("open_inout:mean_history_gatt");
forrest@0	923	exit(2);
forrest@0	924	}
forrest@0	925	if (!(stddev_history_gatt = (char )malloc((strlen(stddev_prefix) + strlen(flist)+1)sizeof(char)))) {
forrest@0	926	perror("open_inout:stddev_history_gatt");
forrest@0	927	exit(2);
forrest@0	928	}
forrest@0	929	sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
forrest@0	930	sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
forrest@0	931
forrest@1	932	/* The two time_bounds values must be handled explicitly by obtaining
forrest@1	933	* the first one from the first file and the last one from the last
forrest@1	934	* file. These are then passed to copy_metadata() below. */
forrest@1	935	get_time_bounds(first_fname, last_fname, tbounds);
forrest@1	936	#ifdef DEBUG
forrest@1	937	fprintf(stderr, "Got back tbounds=%lf,%lf\n", tbounds[0], tbounds[1]);
forrest@1	938	#endif /* DEBUG */
forrest@1	939
forrest@1	940	/* Open last input file */
forrest@1	941	wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
forrest@0	942	/* Inquire about number of dimensions, variables, global attributes
forrest@0	943	* and the ID of the unlimited dimension
forrest@0	944	*/
forrest@0	945	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
forrest@0	946	&input_unlimdimid));
forrest@0	947	/* Grab dimension IDs and lengths from input file */
forrest@0	948	get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
forrest@0	949	/* Grab desired global attributes from input file */
forrest@0	950	get_gatts(input_ncid, input_ngatts, &input_gatt_head);
forrest@0	951	/* Grab variable IDs and attributes from input file */
forrest@0	952	get_vars(input_ncid, input_nvars, &input_var_head);
forrest@0	953	#ifdef DEBUG
forrest@0	954	varlist_print(input_var_head);
forrest@0	955	#endif /* DEBUG */
forrest@0	956	/* Create netCDF files for monthly means and standard deviations */
forrest@0	957	/* Create new netCDF files */
forrest@0	958	wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
forrest@0	959	wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
forrest@0	960	/* Define dimensions */
forrest@0	961	mean_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0	962	input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
forrest@0	963	&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
forrest@0	964	stddev_ndims = put_dimensions(input_dim_head, input_ndims,
forrest@0	965	input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
forrest@0	966	&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
forrest@0	967	/* Define variables and their attributes */
forrest@0	968	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0	969	idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
forrest@0	970	MEAN_TYPE);
forrest@0	971	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
forrest@0	972	idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
forrest@0	973	STDDEV_TYPE);
forrest@0	974	/* Store global attributes */
forrest@0	975	mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
forrest@0	976	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
forrest@0	977	stddev_history_gatt);
forrest@0	978	/* End define mode */
forrest@0	979	wrap_nc(nc_enddef(mean_ncid));
forrest@0	980	wrap_nc(nc_enddef(stddev_ncid));
forrest@0	981	/* Write out metdata variables that do not depend on "time" */
forrest@0	982	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
forrest@1	983	mean_var_head, tbounds);
forrest@0	984	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
forrest@1	985	stddev_var_head, tbounds);
forrest@0	986
forrest@0	987	wrap_nc(nc_close(input_ncid));
forrest@0	988
forrest@0	989	return;
forrest@0	990	}
forrest@0	991
forrest@0	992	size_t count_samples(int nifnames, char **input_fnames)
forrest@0	993	{
forrest@0	994	int i;
forrest@0	995	char name[NC_MAX_NAME+1];
forrest@0	996	size_t len, cnt = 0;
forrest@0	997
forrest@0	998	for (i = 0; i < nifnames; i++) {
forrest@0	999	/* Open input file */
forrest@0	1000	wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0	1001	/*
forrest@0	1002	* Inquire about number of dimensions, variables, global
forrest@0	1003	* attributes and the ID of the unlimited dimension
forrest@0	1004	*/
forrest@0	1005	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0	1006	&input_ngatts, &input_unlimdimid));
forrest@0	1007	wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0	1008	if (strcmp(name, time_name)) {
forrest@0	1009	fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0	1010	exit(4);
forrest@0	1011	}
forrest@0	1012	#ifdef DEBUG
forrest@0	1013	printf("%ld samples in %s\n", (long)len, input_fnames[i]);
forrest@0	1014	#endif /* DEBUG */
forrest@0	1015	wrap_nc(nc_close(input_ncid));
forrest@0	1016	cnt += len;
forrest@0	1017	}
forrest@0	1018	return cnt;
forrest@0	1019	}
forrest@0	1020
forrest@0	1021	void alloc_means_stddevs(size_t d1, size_t d2, double *meansp, double stddevsp, size_t **cell_samplesp)
forrest@0	1022	{
forrest@0	1023	/*
forrest@0	1024	* Allocate space for arrays of means and standard deviations by
forrest@0	1025	* hour of day.
forrest@0	1026	*/
forrest@0	1027	int i;
forrest@0	1028	size_t **cell_samples;
forrest@0	1029	double means, stddevs;
forrest@0	1030
forrest@0	1031	if (!(meansp = (double )malloc(HOURS_PER_DAY sizeof(double *)))) {
forrest@0	1032	perror("alloc_means_stddevs:*meansp");
forrest@0	1033	exit(2);
forrest@0	1034	}
forrest@0	1035	if (!(stddevsp = (double )malloc(HOURS_PER_DAY sizeof(double *)))) {
forrest@0	1036	perror("alloc_means_stddevs:*stddevsp");
forrest@0	1037	exit(2);
forrest@0	1038	}
forrest@0	1039	if (!(cell_samplesp = (size_t )malloc(HOURS_PER_DAY sizeof(size_t *)))) {
forrest@0	1040	perror("alloc_means_stddevs:*cell_samplesp");
forrest@0	1041	exit(2);
forrest@0	1042	}
forrest@0	1043
forrest@0	1044	means = *meansp;
forrest@0	1045	stddevs = *stddevsp;
forrest@0	1046	cell_samples = *cell_samplesp;
forrest@0	1047
forrest@0	1048	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0	1049	if (!(means[i] = (double )malloc(d1 d2 * sizeof(double)))) {
forrest@0	1050	perror("alloc_means_stddevs:means[i]");
forrest@0	1051	exit(2);
forrest@0	1052	}
forrest@0	1053	if (!(stddevs[i] = (double )malloc(d1 d2 * sizeof(double)))) {
forrest@0	1054	perror("alloc_means_stddevs:stddevs[i]");
forrest@0	1055	exit(2);
forrest@0	1056	}
forrest@0	1057	if (!(cell_samples[i] = (size_t )malloc(d1 d2 * sizeof(size_t)))) {
forrest@0	1058	perror("alloc_means_stddevs:cell_samples[i]");
forrest@0	1059	exit(2);
forrest@0	1060	}
forrest@0	1061	}
forrest@0	1062
forrest@0	1063	return;
forrest@0	1064	}
forrest@0	1065
forrest@0	1066	void init_means_stddevs(size_t d1, size_t d2, double means, double stddevs, size_t **cell_samples, float FillValue)
forrest@0	1067	{
forrest@0	1068	int hr, i, j, pos;
forrest@0	1069
forrest@0	1070	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1071	for (i = 0; i < d1; i++) {
forrest@0	1072	#pragma _CRI concurrent
forrest@0	1073	for (j = 0; j < d2; j++) {
forrest@0	1074	pos = i*d2+j;
forrest@0	1075	means[hr][pos] = FillValue;
forrest@0	1076	stddevs[hr][pos] = FillValue;
forrest@0	1077	cell_samples[hr][pos] = 0;
forrest@0	1078	}
forrest@0	1079	}
forrest@0	1080	}
forrest@0	1081
forrest@0	1082	return;
forrest@0	1083	}
forrest@0	1084
forrest@0	1085	void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
forrest@0	1086	{
forrest@0	1087	int i, j, hr, pos;
forrest@0	1088
forrest@0	1089	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1090	for (i = 0; i < d1; i++) {
forrest@0	1091	#pragma _CRI concurrent
forrest@0	1092	for (j = 0; j < d2; j++) {
forrest@0	1093	pos = i*d2+j;
forrest@0	1094	cell_samples[hr][pos] = 0;
forrest@0	1095	}
forrest@0	1096	}
forrest@0	1097	}
forrest@0	1098
forrest@0	1099	return;
forrest@0	1100	}
forrest@0	1101
forrest@0	1102	void free_means_stddevs(double means, double stddevs, size_t **cell_samples)
forrest@0	1103	{
forrest@0	1104	/*
forrest@0	1105	* Free space from arrays of means and standard deviations by
forrest@0	1106	* hour of day.
forrest@0	1107	*/
forrest@0	1108	int i;
forrest@0	1109
forrest@0	1110	for (i = 0; i < HOURS_PER_DAY; i++) {
forrest@0	1111	free(means[i]);
forrest@0	1112	free(stddevs[i]);
forrest@0	1113	free(cell_samples[i]);
forrest@0	1114	}
forrest@0	1115
forrest@0	1116	free(means);
forrest@0	1117	free(stddevs);
forrest@0	1118	free(cell_samples);
forrest@0	1119
forrest@0	1120	return;
forrest@0	1121	}
forrest@0	1122
forrest@0	1123	void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
forrest@0	1124	char FillFlag, float FillValue, double means, size_t cell_samples)
forrest@0	1125	{
forrest@0	1126	int i, j, hr, pos;
forrest@0	1127
forrest@0	1128	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1129
forrest@0	1130	for (i = 0; i < d1; i++) {
forrest@0	1131	#pragma _CRI concurrent
forrest@0	1132	for (j = 0; j < d2; j++) {
forrest@0	1133	pos = i*d2+j;
forrest@0	1134	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue)) {
forrest@0	1135	if (cell_samples[hr][pos] == 0)
forrest@0	1136	means[hr][pos] = (double)val[pos];
forrest@0	1137	else
forrest@0	1138	means[hr][pos] += (double)val[pos];
forrest@0	1139	++cell_samples[hr][pos];
forrest@0	1140	}
forrest@0	1141	}
forrest@0	1142	}
forrest@0	1143
forrest@0	1144	return;
forrest@0	1145	}
forrest@0	1146
forrest@0	1147	void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
forrest@0	1148	char FillFlag, float FillValue, double means, size_t cell_samples)
forrest@0	1149	{
forrest@0	1150	int i, j, hr, pos;
forrest@0	1151
forrest@0	1152	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1153
forrest@0	1154	for (i = 0; i < d1; i++) {
forrest@0	1155	#pragma _CRI concurrent
forrest@0	1156	for (j = 0; j < d2; j++) {
forrest@0	1157	pos = i*d2+j;
forrest@0	1158	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue)) {
forrest@0	1159	if (cell_samples[hr][pos] == 0)
forrest@0	1160	means[hr][pos] = val[pos];
forrest@0	1161	else
forrest@0	1162	means[hr][pos] += val[pos];
forrest@0	1163	++cell_samples[hr][pos];
forrest@0	1164	}
forrest@0	1165	}
forrest@0	1166	}
forrest@0	1167
forrest@0	1168	return;
forrest@0	1169	}
forrest@0	1170
forrest@0	1171
forrest@0	1172	void divide_means(size_t d1, size_t d2, double means, size_t cell_samples)
forrest@0	1173	{
forrest@0	1174	int i, j, hr, pos;
forrest@0	1175
forrest@0	1176	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1177	for (i = 0; i < d1; i++) {
forrest@0	1178	#pragma _CRI concurrent
forrest@0	1179	for (j = 0; j < d2; j++) {
forrest@0	1180	pos = i*d2+j;
forrest@0	1181	if (cell_samples[hr][pos] != 0) {
forrest@0	1182	means[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0	1183	}
forrest@0	1184	}
forrest@0	1185	}
forrest@0	1186	}
forrest@0	1187
forrest@0	1188	return;
forrest@0	1189	}
forrest@0	1190
forrest@0	1191	void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
forrest@0	1192	char FillFlag, float FillValue, double means, double stddevs,
forrest@0	1193	size_t **cell_samples)
forrest@0	1194	{
forrest@0	1195	int i, j, hr, pos;
forrest@0	1196	double delta;
forrest@0	1197
forrest@0	1198	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1199
forrest@0	1200	for (i = 0; i < d1; i++) {
forrest@0	1201	#pragma _CRI concurrent
forrest@0	1202	for (j = 0; j < d2; j++) {
forrest@0	1203	pos = i*d2+j;
forrest@0	1204	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue
forrest@0	1205	&& means[hr][pos] != FillValue)) {
forrest@0	1206	delta = means[hr][pos] - (double)val[pos];
forrest@0	1207	if (cell_samples[hr][pos] == 0)
forrest@0	1208	stddevs[hr][pos] = delta * delta;
forrest@0	1209	else
forrest@0	1210	stddevs[hr][pos] += delta * delta;
forrest@0	1211	++cell_samples[hr][pos];
forrest@0	1212	}
forrest@0	1213	}
forrest@0	1214	}
forrest@0	1215
forrest@0	1216	return;
forrest@0	1217	}
forrest@0	1218
forrest@0	1219	void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
forrest@0	1220	char FillFlag, float FillValue, double means, double stddevs,
forrest@0	1221	size_t **cell_samples)
forrest@0	1222	{
forrest@0	1223	int i, j, hr, pos;
forrest@0	1224	double delta;
forrest@0	1225
forrest@0	1226	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
forrest@0	1227
forrest@0	1228	for (i = 0; i < d1; i++) {
forrest@0	1229	#pragma _CRI concurrent
forrest@0	1230	for (j = 0; j < d2; j++) {
forrest@0	1231	pos = i*d2+j;
forrest@0	1232	if (!FillFlag \|\| (FillFlag && val[pos] != FillValue
forrest@0	1233	&& means[hr][pos] != FillValue)) {
forrest@0	1234	delta = means[hr][pos] - val[pos];
forrest@0	1235	if (cell_samples[hr][pos] == 0)
forrest@0	1236	stddevs[hr][pos] = delta * delta;
forrest@0	1237	else
forrest@0	1238	stddevs[hr][pos] += delta * delta;
forrest@0	1239	++cell_samples[hr][pos];
forrest@0	1240	}
forrest@0	1241	}
forrest@0	1242	}
forrest@0	1243
forrest@0	1244	return;
forrest@0	1245	}
forrest@0	1246
forrest@0	1247
forrest@0	1248	void divide_sqrt_stddevs(size_t d1, size_t d2, double stddevs, size_t cell_samples)
forrest@0	1249	{
forrest@0	1250	int i, j, hr, pos;
forrest@0	1251
forrest@0	1252	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1253	for (i = 0; i < d1; i++) {
forrest@0	1254	#pragma _CRI concurrent
forrest@0	1255	for (j = 0; j < d2; j++) {
forrest@0	1256	pos = i*d2+j;
forrest@0	1257	if (cell_samples[hr][pos] != 0) {
forrest@0	1258	stddevs[hr][pos] /= (double)cell_samples[hr][pos];
forrest@0	1259	stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
forrest@0	1260	}
forrest@0	1261	}
forrest@0	1262	}
forrest@0	1263	}
forrest@0	1264
forrest@0	1265	return;
forrest@0	1266	}
forrest@0	1267
forrest@0	1268	void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void *valsp, int mcsecp)
forrest@0	1269	{
forrest@0	1270	int i;
forrest@0	1271	void **vals;
forrest@0	1272
forrest@0	1273	/* Allocate space for all timeslices */
forrest@0	1274	if (!(valsp = (void )malloc(nsamples sizeof(void *)))) {
forrest@0	1275	perror("alloc_all_samples:*valsp");
forrest@0	1276	exit(2);
forrest@0	1277	}
forrest@0	1278	vals = *valsp;
forrest@0	1279	for (i = 0; i < nsamples; i++) {
forrest@0	1280	vals[i] = alloc_var(nctype, (d1 * d2));
forrest@0	1281	}
forrest@0	1282	if (!(mcsecp = (int )malloc(nsamples * sizeof(int)))) {
forrest@0	1283	perror("alloc_all_samples:*mcsecp");
forrest@0	1284	exit(2);
forrest@0	1285	}
forrest@0	1286
forrest@0	1287	return;
forrest@0	1288	}
forrest@0	1289
forrest@0	1290	void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void *vals, int mcsec)
forrest@0	1291	{
forrest@0	1292	int i, j;
forrest@0	1293	float **fvals = NULL;
forrest@0	1294	double **dvals = NULL;
forrest@0	1295
forrest@0	1296	switch(nctype) {
forrest@0	1297	case NC_FLOAT:
forrest@0	1298	fvals = (float **)vals;
forrest@0	1299	break;
forrest@0	1300	case NC_DOUBLE:
forrest@0	1301	dvals = (double **)vals;
forrest@0	1302	break;
forrest@0	1303	default:
forrest@0	1304	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1305	exit(3);
forrest@0	1306	}
forrest@0	1307
forrest@0	1308	for (i = 0; i < nsamples; i++) {
forrest@0	1309	for (j = 0; j < (d1 * d2); j++) {
forrest@0	1310	switch(nctype) {
forrest@0	1311	case NC_FLOAT:
forrest@0	1312	fvals[i][j] = (float)0.0;
forrest@0	1313	break;
forrest@0	1314	case NC_DOUBLE:
forrest@0	1315	dvals[i][j] = (double)0.0;
forrest@0	1316	break;
forrest@0	1317	default:
forrest@0	1318	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1319	exit(3);
forrest@0	1320	}
forrest@0	1321	}
forrest@0	1322	mcsec[i] = 0;
forrest@0	1323	}
forrest@0	1324
forrest@0	1325	return;
forrest@0	1326	}
forrest@0	1327
forrest@0	1328	void free_all_samples(size_t nsamples, void *vals, int mcsec)
forrest@0	1329	{
forrest@0	1330	int i;
forrest@0	1331
forrest@0	1332	for (i = 0; i < nsamples; i++)
forrest@0	1333	if (vals[i]) free(vals[i]);
forrest@0	1334
forrest@0	1335	free(vals);
forrest@0	1336	free(mcsec);
forrest@0	1337
forrest@0	1338	return;
forrest@0	1339	}
forrest@0	1340
forrest@0	1341	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	1342	{
forrest@0	1343	int i, ts;
forrest@0	1344	int varid;
forrest@0	1345	size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	1346	char name[NC_MAX_NAME+1];
forrest@0	1347
forrest@0	1348	for (i = 0; i < nifnames; i++) {
forrest@0	1349	#ifdef DEBUG
forrest@0	1350	printf("\tOpening %s", input_fnames[i]);
forrest@0	1351	if (ndims > 3) printf(" to retrieve level %d\n", level);
forrest@0	1352	#endif /* DEBUG */
forrest@0	1353	/* Open input file */
forrest@0	1354	wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
forrest@0	1355	/*
forrest@0	1356	* Inquire about number of dimensions, variables, global
forrest@0	1357	* attributes and the ID of the unlimited dimension
forrest@0	1358	*/
forrest@0	1359	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
forrest@0	1360	&input_ngatts, &input_unlimdimid));
forrest@0	1361	wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
forrest@0	1362	if (strcmp(name, time_name)) {
forrest@0	1363	fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
forrest@0	1364	exit(4);
forrest@0	1365	}
forrest@0	1366	/* Make sure we don't run off the end of allocated space */
forrest@0	1367	if ((nslice+len) > nsamples) {
forrest@0	1368	fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);
forrest@0	1369	exit(4);
forrest@0	1370	}
forrest@0	1371	/* Get seconds of day variable */
forrest@0	1372	wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
forrest@0	1373	wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));
forrest@0	1374	/* Get variable ID for requested variable */
forrest@0	1375	wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
forrest@0	1376	/* Retrieve time slice of desired variable */
forrest@0	1377	for (ts = 0; ts < len; ts++) {
forrest@0	1378	start[0] = ts;
forrest@0	1379	count[0] = 1;
forrest@0	1380	start[(ndims-2)] = 0;
forrest@0	1381	count[(ndims-2)] = d1;
forrest@0	1382	start[(ndims-1)] = 0;
forrest@0	1383	count[(ndims-1)] = d2;
forrest@0	1384	if (ndims == 4) {
forrest@0	1385	start[1] = level;
forrest@0	1386	count[1] = 1;
forrest@0	1387	}
forrest@0	1388	switch(nctype) {
forrest@0	1389	case NC_FLOAT:
forrest@0	1390	wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));
forrest@0	1391	break;
forrest@0	1392	case NC_DOUBLE:
forrest@0	1393	wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));
forrest@0	1394	break;
forrest@0	1395	default:
forrest@0	1396	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1397	exit(3);
forrest@0	1398	}
forrest@0	1399	nslice++;
forrest@0	1400	}
forrest@0	1401
forrest@0	1402	/* Close input file */
forrest@0	1403	wrap_nc(nc_close(input_ncid));
forrest@0	1404	}
forrest@0	1405
forrest@0	1406	return;
forrest@0	1407	}
forrest@0	1408
forrest@0	1409	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	1410	{
forrest@0	1411	int i;
forrest@0	1412
forrest@0	1413	for (i = 0; i < nsamples; i++) {
forrest@0	1414	switch(nctype) {
forrest@0	1415	case NC_FLOAT:
forrest@0	1416	if (stat_type == MEAN_TYPE)
forrest@0	1417	add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0	1418	else
forrest@0	1419	add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0	1420	break;
forrest@0	1421	case NC_DOUBLE:
forrest@0	1422	if (stat_type == MEAN_TYPE)
forrest@0	1423	add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
forrest@0	1424	else
forrest@0	1425	add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
forrest@0	1426	break;
forrest@0	1427	default:
forrest@0	1428	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1429	exit(3);
forrest@0	1430	}
forrest@0	1431	}
forrest@0	1432
forrest@0	1433	/* Divide sums by number of samples */
forrest@0	1434	if (stat_type == MEAN_TYPE)
forrest@0	1435	divide_means(d1, d2, means, cell_samples);
forrest@0	1436	else
forrest@0	1437	divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
forrest@0	1438
forrest@0	1439	return;
forrest@0	1440	}
forrest@0	1441
forrest@0	1442	float double_to_float(size_t len, double dvar)
forrest@0	1443	{
forrest@0	1444	int i;
forrest@0	1445	float *fvar;
forrest@0	1446
forrest@0	1447	if (!(fvar = (float )malloc(len sizeof(float)))) {
forrest@0	1448	perror("double_to_float:fvar");
forrest@0	1449	exit(2);
forrest@0	1450	}
forrest@0	1451
forrest@0	1452	for (i = 0; i < len; i++)
forrest@0	1453	fvar[i] = (float)dvar[i];
forrest@0	1454
forrest@0	1455	return fvar;
forrest@0	1456	}
forrest@0	1457
forrest@0	1458	void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
forrest@0	1459	int level, int ndims, double **var_hours)
forrest@0	1460	{
forrest@0	1461	/* Output dimensions should be one larger than input dimensions */
forrest@0	1462	int i, hr;
forrest@0	1463	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0	1464	float *fvar = NULL;
forrest@0	1465
forrest@0	1466	if (nctype == NC_FLOAT) {
forrest@0	1467	if (!(fvar = (float )malloc(d1 d2 * sizeof(float)))) {
forrest@0	1468	perror("write_var_hours:fvar");
forrest@0	1469	exit(2);
forrest@0	1470	}
forrest@0	1471	}
forrest@0	1472
forrest@0	1473	for (hr = 0; hr < HOURS_PER_DAY; hr++) {
forrest@0	1474	start[0] = 0;
forrest@0	1475	count[0] = 1; /* time */
forrest@0	1476	start[1] = hr;
forrest@0	1477	count[1] = 1; /* hour */
forrest@0	1478	start[(ndims-2)] = 0;
forrest@0	1479	count[(ndims-2)] = d1;
forrest@0	1480	start[(ndims-1)] = 0;
forrest@0	1481	count[(ndims-1)] = d2;
forrest@0	1482	if (ndims == 5) {
forrest@0	1483	start[2] = level;
forrest@0	1484	count[2] = 1;
forrest@0	1485	}
forrest@0	1486	switch (nctype) {
forrest@0	1487	case NC_FLOAT:
forrest@0	1488	for (i = 0; i < (d1 * d2); i++)
forrest@0	1489	fvar[i] = (float)var_hours[hr][i];
forrest@0	1490	wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
forrest@0	1491	break;
forrest@0	1492	case NC_DOUBLE:
forrest@0	1493	wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
forrest@0	1494	break;
forrest@0	1495	default:
forrest@0	1496	fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
forrest@0	1497	exit(3);
forrest@0	1498	}
forrest@0	1499	}
forrest@0	1500
forrest@0	1501	if (nctype == NC_FLOAT)
forrest@0	1502	free(fvar);
forrest@0	1503
forrest@0	1504	return;
forrest@0	1505	}
forrest@0	1506
forrest@0	1507	void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
forrest@0	1508	{
forrest@0	1509	int j, k, nlevels, *mcsec;
forrest@0	1510	size_t d1, d2, **cell_samples;
forrest@0	1511	double means, stddevs;
forrest@0	1512	struct var in_vnode, out_vnode;
forrest@0	1513	void **vals;
forrest@0	1514
forrest@0	1515	if (input_var_head) {
forrest@0	1516	in_vnode = input_var_head;
forrest@0	1517	/* March through non-metadata variables to compute stats */
forrest@0	1518	for (j = 0; j == 0 \|\| in_vnode != input_var_head; j++) {
forrest@0	1519	if (!in_vnode->metadata) {
forrest@0	1520	/* Make sure time is really the first dimension */
forrest@0	1521	if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
forrest@0	1522	fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
forrest@0	1523	exit(3);
forrest@0	1524	}
forrest@0	1525	/* Figure out input dimensions */
forrest@0	1526	if (in_vnode->ndims < 3 \|\| in_vnode->ndims > 4) {
forrest@0	1527	fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
forrest@0	1528	exit(3);
forrest@0	1529	}
forrest@0	1530	/* Assume 2D output; find dimensions */
forrest@0	1531	d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
forrest@0	1532	d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
forrest@0	1533	if (in_vnode->ndims == 3)
forrest@0	1534	nlevels = 1;
forrest@0	1535	else
forrest@0	1536	nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
forrest@0	1537	/* Allocate working space for means and stddevs */
forrest@0	1538	alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
forrest@0	1539	init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
forrest@0	1540	/* Allocate and initize space for entire field across all files */
forrest@0	1541	alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);
forrest@0	1542	init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);
forrest@0	1543	printf("Computing statistics for %s\n",
forrest@0	1544	in_vnode->name);
forrest@0	1545	/* Compute means and stddevs, then write them */
forrest@0	1546	for (k = 0; k < nlevels; k++) {
forrest@0	1547	/* Read the entire fields from all files */
forrest@0	1548	read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);
forrest@0	1549	/* Compute means */
forrest@0	1550	compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
forrest@0	1551	/* Find corresponding output variable on the mean output file */
forrest@0	1552	out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
forrest@0	1553	/* Write out the means for this variable */
forrest@0	1554	write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
forrest@0	1555	/* Zero out cell_samples so they can be used as a flag again for computing stddevs */
forrest@0	1556	reset_cell_samples(d1, d2, cell_samples);
forrest@0	1557	/* Compute stddevs using means */
forrest@0	1558	compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
forrest@0	1559	/* Find corresponding output variable on the stddev output file */
forrest@0	1560	out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
forrest@0	1561	/* Write out stddevs for this variable */
forrest@0	1562	write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
forrest@0	1563	}
forrest@0	1564
forrest@0	1565	/* Free all samples */
forrest@0	1566	free_all_samples(nsamples, vals, mcsec);
forrest@0	1567	/* Free working space for means and stddevs */
forrest@0	1568	free_means_stddevs(means, stddevs, cell_samples);
forrest@0	1569	}
forrest@0	1570	in_vnode = in_vnode->next;
forrest@0	1571	}
forrest@0	1572	}
forrest@0	1573	return;
forrest@0	1574	}
forrest@0	1575
forrest@0	1576	void usage(char *program)
forrest@0	1577	{
forrest@0	1578	fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
forrest@0	1579	fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
forrest@0	1580	exit(4);
forrest@0	1581	}
forrest@0	1582
forrest@0	1583	int main(int argc, char **argv)
forrest@0	1584	{
forrest@0	1585	int i, ic, nifnames;
forrest@0	1586	size_t len, pos, nsamples;
forrest@0	1587	char mfname, sfname, *ifnames, flist;
forrest@0	1588
forrest@0	1589	ifnames = NULL;
forrest@0	1590	mfname = sfname = flist = NULL;
forrest@0	1591	nifnames = 0;
forrest@0	1592
forrest@0	1593	#ifdef DEBUG
forrest@0	1594	printf("Number of metadata variables (nmvars) = %d\n", nmvars);
forrest@0	1595	fflush(stdout);
forrest@0	1596	#endif /* DEBUG */
forrest@0	1597
forrest@0	1598	/* check command-line flags and switches */
forrest@0	1599	while ((ic = getopt(argc, argv, "m:s:")) != -1) {
forrest@0	1600	switch(ic) {
forrest@0	1601	case 'm':
forrest@0	1602	if (!(mfname = strdup(optarg))) {
forrest@0	1603	perror("mfname");
forrest@0	1604	exit(2);
forrest@0	1605	}
forrest@0	1606	break;
forrest@0	1607	case 's':
forrest@0	1608	if (!(sfname = strdup(optarg))) {
forrest@0	1609	perror("sfname");
forrest@0	1610	exit(2);
forrest@0	1611	}
forrest@0	1612	break;
forrest@0	1613	}
forrest@0	1614	}
forrest@0	1615
forrest@0	1616	if (!mfname) {
forrest@0	1617	fprintf(stderr, "Output file name for writing means required.\n");
forrest@0	1618	usage(argv[0]);
forrest@0	1619	}
forrest@0	1620	if (!sfname) {
forrest@0	1621	fprintf(stderr, "Output file name for writing standard deviations required.\n");
forrest@0	1622	usage(argv[0]);
forrest@0	1623	}
forrest@0	1624
forrest@0	1625	if (optind < argc) {
forrest@0	1626	if (!(ifnames = (char *)malloc((argc-optind+1)sizeof(char *)))) {
forrest@0	1627	perror("ifnames");
forrest@0	1628	exit(2);
forrest@0	1629	}
forrest@0	1630	for (i = optind; i < argc; i++) {
forrest@0	1631	if (!(ifnames[nifnames++] = strdup(argv[i]))) {
forrest@0	1632	perror("ifnames[nifnames++]");
forrest@0	1633	exit(2);
forrest@0	1634	}
forrest@0	1635	}
forrest@0	1636	}
forrest@0	1637	else {
forrest@0	1638	fprintf(stderr, "At least one input file name is required.\n");
forrest@0	1639	usage(argv[0]);
forrest@0	1640	}
forrest@0	1641
forrest@0	1642
forrest@0	1643	/*
forrest@0	1644	* Build list of input files to be included in the global history
forrest@0	1645	* attribute on the two outputs files.
forrest@0	1646	*/
forrest@0	1647	for (i = len = 0; i < nifnames; i++)
forrest@0	1648	len += strlen(ifnames[i]);
forrest@0	1649	len += nifnames + 1; /* space in front of every name + null terminator */
forrest@0	1650	if (!(flist = (char )malloc(len sizeof(char)))) {
forrest@0	1651	perror("flist");
forrest@0	1652	exit(2);
forrest@0	1653	}
forrest@0	1654	for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
forrest@0	1655	sprintf(&flist[pos], " %s", ifnames[i]);
forrest@0	1656	#ifdef DEBUG
forrest@0	1657	printf("flist=%s\n", flist);
forrest@0	1658	fflush(stdout);
forrest@0	1659	#endif /* DEBUG */
forrest@0	1660
forrest@0	1661	/* Open every file to count up the number of time samples. */
forrest@0	1662	nsamples = count_samples(nifnames, ifnames);
forrest@0	1663	#ifdef DEBUG
forrest@0	1664	printf("Number of samples across all files = %ld\n", (long)nsamples);
forrest@0	1665	fflush(stdout);
forrest@0	1666	#endif /* DEBUG */
forrest@0	1667
forrest@0	1668	/*
forrest@0	1669	* Open the last input file and the two output files (for mean and
forrest@0	1670	* standard deviation). Define dimensions, variables, and attributes
forrest@0	1671	* for the two output files based on the last input files. The
forrest@0	1672	* last file is used because some metadata variables will contain
forrest@0	1673	* only values from the last time slice from the period over which
forrest@0	1674	* the statistics are computed.
forrest@0	1675	*/
forrest@1	1676	open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
forrest@1	1677	nsamples);
forrest@0	1678
forrest@0	1679	compute_stats(nifnames, ifnames, nsamples);
forrest@0	1680
forrest@0	1681	/* Close the two output files */
forrest@0	1682	wrap_nc(nc_close(mean_ncid));
forrest@0	1683	wrap_nc(nc_close(stddev_ncid));
forrest@0	1684
forrest@0	1685	return 0;
forrest@0	1686	}

author	Forrest Hoffman <forrest@climatemodeling.org>
	Wed, 10 Oct 2007 11:59:02 -0400
changeset 4	dd8e6719647b
parent 1	2ce4ee911439
permissions	-rw-r--r--