c-lamp_postproc: h1_summary2.c@dd8e6719647b

     1 #include <stdio.h>

     2 #include <unistd.h>

     3 #include <stdlib.h>

     4 #include <string.h>

     5 #include <math.h>

     6 #include "netcdf.h"

     8 /*

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

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

    11  * by hour of day.

    12  */

    14 #define	MEAN_TYPE	0

    15 #define STDDEV_TYPE	1

    17 #define HOURS_PER_DAY	24

    18 #define MIN_PER_HOUR	60

    19 #define SEC_PER_MIN	60

    20 #define SEC_PER_HOUR	(MIN_PER_HOUR * SEC_PER_MIN)

    22 static char *metadata_vars[] = {

    23 	"lon",

    24 	"lat",

    25 	"lonrof",

    26 	"latrof",

    27 	"time",

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

    29 	"levsoi",

    30 	"levlak",

    31 	"edgen",

    32 	"edgee",

    33 	"edges",

    34 	"edgew",

    35 	"longxy",

    36 	"latixy",

    37 	"area",

    38 	"areaupsc",

    39 	"topo",

    40 	"topodnsc",

    41 	"landfrac",

    42 	"landmask",

    43 	"pftmask",

    44 	"indxupsc",

    45 	"mcdate",

    46 	"mcsec",

    47 	"mdcur",

    48 	"mscur",

    49 	"nstep",

    50 	"time_bounds",

    51 	"date_written",

    52 	"time_written"

    53 };

    55 struct text_att {

    56 	char *name;

    57 	char *value;

    58 	struct text_att *next;

    59 	struct text_att *prev;

    60 };

    62 struct dim {

    63 	int dimid;

    64 	char *name;

    65 	size_t len;

    66 	struct dim *next;

    67 	struct dim *prev;

    68 };

    70 struct var {

    71 	int ncvarid;

    72 	char *name;

    73 	nc_type nctype;

    74 	int ndims;

    75 	int *dimids;

    76 	int natts;

    77 	char metadata;

    78 	char FillFlag;

    79 	float FillValue;

    80 	struct var *next;

    81 	struct var *prev;

    82 };

    84 static char *time_name = "time";

    85 static char *time_bounds_name = "time_bounds";

    86 static char *mcsec_name = "mcsec";

    87 static char *history_name = "history";

    88 static char *nsamples_name = "nsamples";

    89 static char *hour_name = "hour", *hour_long_name = "hour of day",

    90 	*hour_units = "hour";

    91 static char *cell_method_name = "cell_method";

    92 /* input stuff */

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

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

    95 static struct text_att *input_gatt_head = NULL;

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

    97 static struct var *input_var_head = NULL;

    98 /* translation stuff */

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

   100 /* output stuff */

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

   102 static int mean_hour_dimid; /* special notes */

   103 static struct dim *mean_dim_head = NULL;

   104 static struct var *mean_var_head = NULL;

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

   106 static int stddev_hour_dimid; /* special notes */

   107 static struct dim *stddev_dim_head = NULL;

   108 static struct var *stddev_var_head = NULL;

   110 char is_metadata(char *name)

   111 {

   112 	int i;

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

   116 	if (i < nmvars)

   117 		return 1;

   118 	else

   119 		return 0;

   120 }

   122 #if 0

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

   124 {

   125 	int i;

   126 	struct dim *p = NULL;

   128 	if (head) {

   129 		node = head;

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

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

   132 				p = node;

   133 			node = node->next;

   134 		}

   135 		return p;

   136 	}

   137 	else

   138 		return NULL;

   139 }

   140 #endif

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

   143 {

   144 	int i;

   145 	struct var *node;

   147 	if (head) {

   148 		node = head;

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

   150 			node = node->next;

   151 		if (i && node == head)

   152 			return NULL;

   153 		else

   154 			return node;

   155 	}

   156 	else

   157 		return NULL;

   158 }

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

   161 {

   162 	struct text_att *head, *node;

   164 	head = *headp;

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

   167 		perror("gattlist_add_node");

   168 		exit(2);

   169 	}

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

   172 		perror("gattlist_add_node");

   173 		exit(2);

   174 	}

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

   176 		perror("gattlist_add_node");

   177 		exit(2);

   178 	}

   180 	if (head == NULL) {

   181 		/* first one */

   182 		*headp = node;

   183 		node->next = node;

   184 		node->prev = node;

   185 	}

   186 	else {

   187 		node->next = head;

   188 		node->prev = head->prev;

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

   190 		head->prev = node;

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

   192 		node->prev->next = node;

   193 	}

   195 	return;

   196 }

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

   199 {

   200 	struct dim *head, *node;

   202 	head = *headp;

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

   205 		perror("dimlist_add_node");

   206 		exit(2);

   207 	}

   209 	node->dimid = dimid;

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

   211 		perror("dimlist_add_node");

   212 		exit(2);

   213 	}

   214 	node->len = len;

   216 	if (head == NULL) {

   217 		/* first one */

   218 		*headp = node;

   219 		node->next = node;

   220 		node->prev = node;

   221 	}

   222 	else {

   223 		node->next = head;

   224 		node->prev = head->prev;

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

   226 		head->prev = node;

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

   228 		node->prev->next = node;

   229 	}

   231 	return node;

   232 }

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

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

   236 	float FillValue)

   237 {

   238 	int i;

   239 	struct var *head, *node;

   241 	head = *headp;

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

   244 		perror("varlist_add_node");

   245 		exit(3);

   246 	}

   248 	node->ncvarid = ncvarid;

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

   250 		perror("varlist_add_node");

   251 		exit(3);

   252 	}

   253 	node->nctype = nctype;

   254 	node->ndims = ndims;

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

   256 		perror("varlist_add_node");

   257 		exit(3);

   258 	}

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

   260 	node->natts = natts;

   261 	node->metadata = is_metadata(name);

   262 	node->FillFlag = FillFlag;

   263 	node->FillValue = FillValue;

   265 	if (head == NULL) {

   266 		/* first one */

   267 		*headp = node;

   268 		node->next = node;

   269 		node->prev = node;

   270 	}

   271 	else {

   272 		node->next = head;

   273 		node->prev = head->prev;

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

   275 		head->prev = node;

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

   277 		node->prev->next = node;

   278 	}

   280 	return;

   281 }

   283 void varlist_print(struct var *headp)

   284 {

   285 	int i, j;

   286 	struct var *node;

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

   290 	if (headp) {

   291 		node = headp;

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

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

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

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

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

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

   298 			printf("\tdimids =");

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

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

   301 			printf("\n");

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

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

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

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

   306 			node = node->next;

   307 		}

   308 	}

   309 	else {

   310 		printf("\tList undefined\n");

   311 	}

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

   315 	return;

   316 }

   318 void wrap_nc(int status)

   319 {

   320 	if (status != NC_NOERR) {

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

   322 		exit(1);

   323 	}

   325 	return;

   326 }

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

   329 {

   330 	int i;

   331 	char name[NC_MAX_NAME+1];

   332 	size_t len;

   333 	struct dim **in_dim_idx;

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

   336 		perror("*in_dim_idxp");

   337 		exit(3);

   338 	}

   339 	in_dim_idx = *in_dim_idxp;

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

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

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

   344 	}

   346 	return;

   347 }

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

   349 {

   350 	int i;

   351 	char name[NC_MAX_NAME+1], *value;

   352 	nc_type xtype;

   353 	size_t len;

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

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

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

   358 		if (xtype != NC_CHAR) {

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

   360 			exit(2);

   361 		}

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

   363 			perror("get_gatts: value");

   364 			exit(3);

   365 		}

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

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

   368 		gattlist_add_node(gatt_headp, name, value);

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

   370 	}

   372 	return;

   373 }

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

   376 {

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

   378 	size_t len;

   379 	float FillValue;

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

   381 	nc_type xtype, atype;

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

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

   385 			&natts));

   386 		/* Check for _FillValue */

   387 		FillFlag = 0;

   388 		FillValue = 0.;

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

   390 			== NC_NOERR) {

   391 			if (atype == NC_FLOAT && len) {

   392 				wrap_nc(nc_get_att_float(ncid, i,

   393 					fill_att_name, &FillValue));

   394 				FillFlag = 1;

   395 			}

   396 		}

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

   398 			natts, FillFlag, FillValue);

   399 	}

   401 	return;

   402 }

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

   405 	size_t nsamples, int **in2outp, int out_ncid,

   406 	struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)

   407 {

   408 	/*

   409 	 * Define dimensions on new files and build translation tables between

   410 	 * dimension IDs.

   411 	 */

   412 	int j, dimid, ndims, *in2out;

   413 	size_t len;

   414 	struct dim *dnode;

   416 	ndims = 0;

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

   419 		perror("put_dimensions: in2outp");

   420 		exit(3);

   421 	}

   422 	in2out = *in2outp;

   424 	if (in_dim_head) {

   425 		dnode = in_dim_head;

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

   427 			if (dnode->dimid == in_unlimdimid)

   428 				len = NC_UNLIMITED;

   429 			else

   430 				len = dnode->len;

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

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

   433 			++ndims;

   434 			in2out[dnode->dimid] = dimid;

   435 			if (dnode->dimid == in_unlimdimid)

   436 				*out_unlimdimidp = dimid;

   437 			/*

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

   439 			 * "nsamples" and "hour" dimensions.

   440 			 */

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

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

   443 				dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);

   444 				++ndims;

   446 				wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));

   447 				dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);

   448 				++ndims;

   449 				/* Track hour dimid for out files */

   450 				*out_hour_dimidp = dimid;

   451 			}

   453 			dnode = dnode->next;

   454 		}

   455 	}

   456 	else {

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

   458 	}

   460 	return ndims;

   461 }

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

   464 {

   465 	/*

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

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

   468 	 */

   469 	int j, hflag = 0, ngatts;

   470 	char *value;

   471 	struct text_att *anode;

   473 	ngatts = 0;

   475 	if (in_gatt_head) {

   476 		anode = in_gatt_head;

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

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

   479 				value = out_history;

   480 				++hflag;

   481 			}

   482 			else

   483 				value = anode->value;

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

   485 			++ngatts;

   486 			anode = anode->next;

   487 		}

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

   489 		if (!hflag) {

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

   491 			++ngatts;

   492 		}

   493 	}

   494 	else {

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

   496 	}

   498 	return ngatts;

   499 }

   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)

   502 {

   503 	/*

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

   505 	 * For normal time-based variables, add a new dimension for hour of

   506 	 * day.  For metadata variables, do not add new dimensions, even if

   507 	 * it is time-based.  Return (possibly new) number of dimensions.

   508 	 */

   509 	int i, ndims;

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

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

   513 		if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)

   514 			&& !metadata) {

   515 			ndims++;

   516 			out_dimids[ndims] = hour_dimid;

   517 		}

   518 	}

   520 	return ndims;

   521 }

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

   524 	int in_varid, int out_ncid, int out_varid)

   525 {

   526 	/*

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

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

   529 	 * the cell_method attribute is set either to "time: mean" or

   530 	 * "time: standard_deviation", depending on the output file type.

   531 	 */

   533 	int i, natts;

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

   535 	char *cell_methods[] = { "time: mean", "time: standard_deviation" };

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

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

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

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

   541 			cmflag = 1;

   542 		}

   543 		else

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

   545 	}

   546 	/*

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

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

   549 	 * on the output file.

   550 	 */

   551 	if (!cmflag && !metadata) {

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

   553 		natts++;

   554 	}

   556 	return natts;

   557 }

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

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

   561 	int stat_type)

   562 {

   563 	/*

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

   565 	 * Return number of variables defined.

   566 	 */

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

   568 	struct var *vnode;

   570 	nvars = 0;

   572 	if (in_var_head) {

   573 		vnode = in_var_head;

   574 		/*

   575 		 * March through input variables creating (mostly) the same

   576 		 * variables on the output file.

   577 		 */

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

   579 			ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);

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

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

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

   583 			++nvars;

   584 			/*

   585 			 * Just after the "time" variable, add the new "hour"

   586 			 * variable for hour of day.

   587 			 */

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

   589 				ndims = 1;

   590 				dimids[0] = hour_dimid;

   591 				wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));

   592 				wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));

   593 				wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));

   594 				varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);

   595 				++nvars;

   596 			}

   598 			vnode = vnode->next;

   599 		}

   600 	}

   601 	else {

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

   603 	}

   605 	return nvars;

   606 }

   608 void *alloc_var(nc_type nctype, size_t len)

   609 {

   610 	void *val;

   612 	switch (nctype) {

   613 		case NC_FLOAT:

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

   615 				perror("alloc_var: val");

   616 				exit(3);

   617 			}

   618 			break;

   619 		case NC_INT:

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

   621 				perror("alloc_var: val");

   622 				exit(3);

   623 			}

   624 			break;

   625 		case NC_DOUBLE:

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

   627 				perror("alloc_var: val");

   628 				exit(3);

   629 			}

   630 			break;

   631 		case NC_CHAR:

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

   633 				perror("alloc_var: val");

   634 				exit(3);

   635 			}

   636 			break;

   637 		default:

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

   639 			exit(3);

   640 	}

   642 	return val;

   643 }

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

   646 {

   647 	int i;

   648 	size_t len = (size_t)1;

   649 	void *val;

   651 	/* Figure out the total size */

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

   654 	val = alloc_var(nctype, len);

   656 	switch (nctype) {

   657 		case NC_FLOAT:

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

   659 			break;

   660 		case NC_INT:

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

   662 			break;

   663 		case NC_DOUBLE:

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

   665 			break;

   666 		case NC_CHAR:

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

   668 			break;

   669 		default:

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

   671 			exit(3);

   672 	}

   674 	return val;

   675 }

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

   678 {

   679 	int i;

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

   681 	void *val;

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

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

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

   686 		exit(3);

   687 	}

   688 	/*

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

   690 	 * and build start/count vectors.

   691 	 */

   692 	start[0] = tslice;

   693 	count[0] = 1;

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

   695 		start[i] = 0;

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

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

   698 	}

   700 	val = alloc_var(nctype, len);

   702 	switch (nctype) {

   703 		case NC_FLOAT:

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

   705 			break;

   706 		case NC_INT:

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

   708 			break;

   709 		case NC_DOUBLE:

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

   711 			break;

   712 		case NC_CHAR:

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

   714 			break;

   715 		default:

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

   717 			exit(3);

   718 	}

   720 	return val;

   721 }

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

   724 {

   725 	switch (nctype) {

   726 		case NC_FLOAT:

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

   728 			break;

   729 		case NC_INT:

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

   731 			break;

   732 		case NC_DOUBLE:

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

   734 			break;

   735 		case NC_CHAR:

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

   737 			break;

   738 		default:

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

   740 			exit(3);

   741 	}

   743 	return;

   744 }

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

   747 {

   748 	int i;

   749 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

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

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

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

   754 		exit(3);

   755 	}

   757 	/* Build start/count vectors */

   758 	start[0] = tslice;

   759 	count[0] = 1;

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

   761 		start[i] = 0;

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

   763 	}

   765 	switch (nctype) {

   766 		case NC_FLOAT:

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

   768 			break;

   769 		case NC_INT:

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

   771 			break;

   772 		case NC_DOUBLE:

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

   774 			break;

   775 		case NC_CHAR:

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

   777 			break;

   778 		default:

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

   780 			exit(3);

   781 	}

   783 	return;

   784 }

   786 void copy_metadata(int in_ncid, struct var *in_var_head,

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

   788 	double *tbounds)

   789 {

   790 	int i, j;

   791 	float hr[HOURS_PER_DAY], mean_time;

   792 	struct var *in_vnode, *out_vnode;

   793 	void *val;

   795 	for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;

   797 	if (in_var_head) {

   798 		in_vnode = in_var_head;

   799 		/*

   800 		 * March through input variables to stuff metadata values into

   801 		 * the new output files. NOTE: Time-based metadata variables

   802 		 * should contain only the last (time-slice) value (from all

   803 		 * input files).

   804 		 */

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

   806 			if (in_vnode->metadata) {

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

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

   809 					/* time is not the first dimension */

   810 #ifdef DEBUG

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

   812 						in_vnode->name);

   813 #endif /* DEBUG */

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

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

   816 				}

   817 				else {

   818 					/* time is the first dimension */

   819 #ifdef DEBUG

   820 					printf("Copying last value of \

   821 time-based metadata variable %s\n", in_vnode->name);

   822 #endif /* DEBUG */

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

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

   825 						write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, tbounds, 0);

   826 					else if (!strcmp(in_vnode->name, time_name)) {

   827 						/* force the timestamp to be the

   828 						 * mean of the time bounds */

   829 						mean_time = (tbounds[0] + tbounds[1]) / 2.0;

   830 						write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, &mean_time, 0);

   831 					}

   832 					else

   833 						write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);

   834 				}

   835 				free(val);

   836 				/*

   837 				 * Just after the "time" variable, write out

   838 				 * the "hour" variable values.

   839 				 */

   840 				if (!strcmp(in_vnode->name, time_name)) {

   841 					out_vnode = varlist_find_by_name(out_var_head, hour_name);

   842 					write_var(out_ncid, out_vnode->ncvarid,

   843 						out_vnode->nctype, hr);

   844 				}

   845 			}

   846 			else {

   847 #ifdef DEBUG

   848 				printf("Skipping variable %s\n",

   849 					in_vnode->name);

   850 #endif /* DEBUG */

   851 			}

   852 			in_vnode = in_vnode->next;

   853 		}

   854 	}

   856 	return;

   857 }

   859 void get_time_bounds(char *first_fname, char *last_fname, double *tbounds)

   860 {

   861 	int time_dimid, time_bounds_varid;

   862 	size_t time_len, time_bounds_index[2];

   863 	double bnd1, bnd2;

   865 	/* Open first input file */

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

   867 	/* Get dimension ID of the time dimension */

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

   869 	/* Get length of time dimension */

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

   871 	/* Get variable ID of time_bounds variable */

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

   873 	/* Set index for reading the first value of time_bounds from the

   874 	 * first timeslice of the first file */

   875 	time_bounds_index[0] = time_bounds_index[1] = 0;

   876 	/* Read the value */

   877 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,

   878 		time_bounds_index, &bnd1));

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

   880 	 * open the second one */

   881 	if (strcmp(first_fname, last_fname)) {

   882 		/* Close the first input file */

   883 		wrap_nc(nc_close(input_ncid));

   884 		/* Open the last input file */

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

   886 		/* Get dimension ID of the time dimension */

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

   888 		/* Get length of time dimension */

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

   890 		/* Get variable ID of time_bounds variable */

   891 		wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,

   892 			&time_bounds_varid));

   893 	}

   894 	/* Set index for reading the second value of time_bounds from the

   895 	 * last timeslice of the last file */

   896 	time_bounds_index[0] = time_len - 1; time_bounds_index[1] = 1;

   897 	/* Read the value */

   898 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,

   899 		time_bounds_index, &bnd2));

   901 	/* Close the last input file */

   902 	wrap_nc(nc_close(input_ncid));

   904 	tbounds[0] = bnd1;

   905 	tbounds[1] = bnd2;

   907 	return;

   908 }

   910 void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)

   911 {

   912 	char *mean_history_gatt, *stddev_history_gatt,

   913 		*mean_prefix="Statistical means from history files:",

   914 		*stddev_prefix="Statistical standard deviations from history files:";

   915 	double tbounds[2];

   917 	/*

   918 	 * Construct strings for history global attributes for the two output

   919 	 * files.

   920 	 */

   921 	if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {

   922 		perror("open_inout:mean_history_gatt");

   923 		exit(2);

   924 	}

   925 	if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {

   926 		perror("open_inout:stddev_history_gatt");

   927 		exit(2);

   928 	}

   929 	sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);

   930 	sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);

   932 	/* The two time_bounds values must be handled explicitly by obtaining

   933 	 * the first one from the first file and the last one from the last

   934 	 * file.  These are then passed to copy_metadata() below. */

   935 	get_time_bounds(first_fname, last_fname, tbounds);

   936 #ifdef DEBUG

   937 	fprintf(stderr, "Got back tbounds=%lf,%lf\n", tbounds[0], tbounds[1]);

   938 #endif /* DEBUG */

   940 	/* Open last input file */

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

   942 	/* Inquire about number of dimensions, variables, global attributes

   943 	 * and the ID of the unlimited dimension

   944 	 */

   945 	wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,

   946 		&input_unlimdimid));

   947 	/* Grab dimension IDs and lengths from input file */

   948 	get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);

   949 	/* Grab desired global attributes from input file */

   950 	get_gatts(input_ncid, input_ngatts, &input_gatt_head);

   951 	/* Grab variable IDs and attributes from input file */

   952 	get_vars(input_ncid, input_nvars, &input_var_head);

   953 #ifdef DEBUG

   954 	varlist_print(input_var_head);

   955 #endif /* DEBUG */

   956 	/* Create netCDF files for monthly means and standard deviations */

   957 	/* Create new netCDF files */

   958 	wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));

   959 	wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));

   960 	/* Define dimensions */

   961 	mean_ndims = put_dimensions(input_dim_head, input_ndims,

   962 		input_unlimdimid, nsamples, &idid2mdid, mean_ncid,

   963 		&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);

   964 	stddev_ndims = put_dimensions(input_dim_head, input_ndims,

   965 		input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,

   966 		&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);

   967 	/* Define variables and their attributes */

   968 	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   969 		idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,

   970 		MEAN_TYPE);

   971 	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   972 		idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,

   973 		STDDEV_TYPE);

   974 	/* Store global attributes */

   975 	mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);

   976 	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,

   977 		stddev_history_gatt);

   978 	/* End define mode */

   979 	wrap_nc(nc_enddef(mean_ncid));

   980 	wrap_nc(nc_enddef(stddev_ncid));

   981 	/* Write out metdata variables that do not depend on "time" */

   982 	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,

   983 		mean_var_head, tbounds);

   984 	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,

   985 		stddev_var_head, tbounds);

   987 	wrap_nc(nc_close(input_ncid));

   989 	return;

   990 }

   992 size_t count_samples(int nifnames, char **input_fnames)

   993 {

   994 	int i;

   995 	char name[NC_MAX_NAME+1];

   996 	size_t len, cnt = 0;

   998 	for (i = 0; i < nifnames; i++) {

   999 		/* Open input file */

  1000 		wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));

  1001 		/*

  1002 		 * Inquire about number of dimensions, variables, global

  1003 		 * attributes and the ID of the unlimited dimension

  1004 		 */

  1005 		wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,

  1006 			&input_ngatts, &input_unlimdimid));

  1007 		wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));

  1008 		if (strcmp(name, time_name)) {

  1009 			fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);

  1010 			exit(4);

  1011 		}

  1012 #ifdef DEBUG

  1013 		printf("%ld samples in %s\n", (long)len, input_fnames[i]);

  1014 #endif /* DEBUG */

  1015 		wrap_nc(nc_close(input_ncid));

  1016 		cnt += len;

  1017 	}

  1018 	return cnt;

  1019 }

  1021 void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)

  1022 {

  1023 	/*

  1024 	 * Allocate space for arrays of means and standard deviations by

  1025 	 * hour of day.

  1026 	 */

  1027 	int i;

  1028 	size_t **cell_samples;

  1029 	double **means, **stddevs;

  1031 	if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {

  1032 		perror("alloc_means_stddevs:*meansp");

  1033 		exit(2);

  1034 	}

  1035 	if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {

  1036 		perror("alloc_means_stddevs:*stddevsp");

  1037 		exit(2);

  1038 	}

  1039 	if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {

  1040 		perror("alloc_means_stddevs:*cell_samplesp");

  1041 		exit(2);

  1042 	}

  1044 	means = *meansp;

  1045 	stddevs = *stddevsp;

  1046 	cell_samples = *cell_samplesp;

  1048 	for (i = 0; i < HOURS_PER_DAY; i++) {

  1049 		if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {

  1050 			perror("alloc_means_stddevs:means[i]");

  1051 			exit(2);

  1052 		}

  1053 		if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {

  1054 			perror("alloc_means_stddevs:stddevs[i]");

  1055 			exit(2);

  1056 		}

  1057 		if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {

  1058 			perror("alloc_means_stddevs:cell_samples[i]");

  1059 			exit(2);

  1060 		}

  1061 	}

  1063 	return;

  1064 }

  1066 void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)

  1067 {

  1068 	int hr, i, j, pos;

  1070 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {

  1071 		for (i = 0; i < d1; i++) {

  1072 #pragma _CRI concurrent

  1073 			for (j = 0; j < d2; j++) {

  1074 				pos = i*d2+j;

  1075 				means[hr][pos] = FillValue;

  1076 				stddevs[hr][pos] = FillValue;

  1077 				cell_samples[hr][pos] = 0;

  1078 			}

  1079 		}

  1080 	}

  1082 	return;

  1083 }

  1085 void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)

  1086 {

  1087 	int i, j, hr, pos;

  1089 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {

  1090 		for (i = 0; i < d1; i++) {

  1091 #pragma _CRI concurrent

  1092 			for (j = 0; j < d2; j++) {

  1093 				pos = i*d2+j;

  1094 				cell_samples[hr][pos] = 0;

  1095 			}

  1096 		}

  1097 	}

  1099 	return;

  1100 }

  1102 void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)

  1103 {

  1104 	/*

  1105 	 * Free space from arrays of means and standard deviations by

  1106 	 * hour of day.

  1107 	 */

  1108 	int i;

  1110 	for (i = 0; i < HOURS_PER_DAY; i++) {

  1111 		free(means[i]);

  1112 		free(stddevs[i]);

  1113 		free(cell_samples[i]);

  1114 	}

  1116 	free(means);

  1117 	free(stddevs);

  1118 	free(cell_samples);

  1120 	return;

  1121 }

  1123 void add_to_means_float(float *val, int sec, size_t d1, size_t d2,

  1124 	char FillFlag, float FillValue, double **means, size_t **cell_samples)

  1125 {

  1126 	int i, j, hr, pos;

  1128 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);

  1130 	for (i = 0; i < d1; i++) {

  1131 #pragma _CRI concurrent

  1132 		for (j = 0; j < d2; j++) {

  1133 			pos = i*d2+j;

  1134 			if (!FillFlag || (FillFlag && val[pos] != FillValue)) {

  1135 				if (cell_samples[hr][pos] == 0)

  1136 					means[hr][pos] = (double)val[pos];

  1137 				else

  1138 					means[hr][pos] += (double)val[pos];

  1139 				++cell_samples[hr][pos];

  1140 			}

  1141 		}

  1142 	}

  1144 	return;

  1145 }

  1147 void add_to_means_double(double *val, int sec, size_t d1, size_t d2,

  1148 	char FillFlag, float FillValue, double **means, size_t **cell_samples)

  1149 {

  1150 	int i, j, hr, pos;

  1152 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);

  1154 	for (i = 0; i < d1; i++) {

  1155 #pragma _CRI concurrent

  1156 		for (j = 0; j < d2; j++) {

  1157 			pos = i*d2+j;

  1158 			if (!FillFlag || (FillFlag && val[pos] != FillValue)) {

  1159 				if (cell_samples[hr][pos] == 0)

  1160 					means[hr][pos] = val[pos];

  1161 				else

  1162 					means[hr][pos] += val[pos];

  1163 				++cell_samples[hr][pos];

  1164 			}

  1165 		}

  1166 	}

  1168 	return;

  1169 }

  1172 void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)

  1173 {

  1174 	int i, j, hr, pos;

  1176 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {

  1177 		for (i = 0; i < d1; i++) {

  1178 #pragma _CRI concurrent

  1179 			for (j = 0; j < d2; j++) {

  1180 				pos = i*d2+j;

  1181 				if (cell_samples[hr][pos] != 0) {

  1182 					means[hr][pos] /= (double)cell_samples[hr][pos];

  1183 				}

  1184 			}

  1185 		}

  1186 	}

  1188 	return;

  1189 }

  1191 void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,

  1192 	char FillFlag, float FillValue, double **means, double **stddevs,

  1193 	size_t **cell_samples)

  1194 {

  1195 	int i, j, hr, pos;

  1196 	double delta;

  1198 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);

  1200 	for (i = 0; i < d1; i++) {

  1201 #pragma _CRI concurrent

  1202 		for (j = 0; j < d2; j++) {

  1203 			pos = i*d2+j;

  1204 			if (!FillFlag || (FillFlag && val[pos] != FillValue

  1205 				&& means[hr][pos] != FillValue)) {

  1206 				delta = means[hr][pos] - (double)val[pos];

  1207 				if (cell_samples[hr][pos] == 0)

  1208 					stddevs[hr][pos] = delta * delta;

  1209 				else

  1210 					stddevs[hr][pos] += delta * delta;

  1211 				++cell_samples[hr][pos];

  1212 			}

  1213 		}

  1214 	}

  1216 	return;

  1217 }

  1219 void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,

  1220 	char FillFlag, float FillValue, double **means, double **stddevs,

  1221 	size_t **cell_samples)

  1222 {

  1223 	int i, j, hr, pos;

  1224 	double delta;

  1226 	hr = (int)floor((double)sec/(double)SEC_PER_HOUR);

  1228 	for (i = 0; i < d1; i++) {

  1229 #pragma _CRI concurrent

  1230 		for (j = 0; j < d2; j++) {

  1231 			pos = i*d2+j;

  1232 			if (!FillFlag || (FillFlag && val[pos] != FillValue

  1233 				&& means[hr][pos] != FillValue)) {

  1234 				delta = means[hr][pos] - val[pos];

  1235 				if (cell_samples[hr][pos] == 0)

  1236 					stddevs[hr][pos] = delta * delta;

  1237 				else

  1238 					stddevs[hr][pos] += delta * delta;

  1239 				++cell_samples[hr][pos];

  1240 			}

  1241 		}

  1242 	}

  1244 	return;

  1245 }

  1248 void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)

  1249 {

  1250 	int i, j, hr, pos;

  1252 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {

  1253 		for (i = 0; i < d1; i++) {

  1254 #pragma _CRI concurrent

  1255 			for (j = 0; j < d2; j++) {

  1256 				pos = i*d2+j;

  1257 				if (cell_samples[hr][pos] != 0) {

  1258 					stddevs[hr][pos] /= (double)cell_samples[hr][pos];

  1259 					stddevs[hr][pos] = sqrt(stddevs[hr][pos]);

  1260 				}

  1261 			}

  1262 		}

  1263 	}

  1265 	return;

  1266 }

  1268 void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp)

  1269 {

  1270 	int i;

  1271 	void **vals;

  1273 	/* Allocate space for all timeslices */

  1274 	if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) {

  1275 		perror("alloc_all_samples:*valsp");

  1276 		exit(2);

  1277 	}

  1278 	vals = *valsp;

  1279 	for (i = 0; i < nsamples; i++) {

  1280 		vals[i] = alloc_var(nctype, (d1 * d2));

  1281 	}

  1282 	if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) {

  1283 		perror("alloc_all_samples:*mcsecp");

  1284 		exit(2);

  1285 	}

  1287 	return;

  1288 }

  1290 void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec)

  1291 {

  1292 	int i, j;

  1293         float **fvals = NULL;

  1294 	double **dvals = NULL;

  1296 	switch(nctype) {

  1297 		case NC_FLOAT:

  1298 			fvals = (float **)vals;

  1299 			break;

  1300 		case NC_DOUBLE:

  1301 			dvals = (double **)vals;

  1302 			break;

  1303 		default:

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

  1305 			exit(3);

  1306 	}

  1308 	for (i = 0; i < nsamples; i++) {

  1309 		for (j = 0; j < (d1 * d2); j++) {

  1310 			switch(nctype) {

  1311 				case NC_FLOAT:

  1312 					fvals[i][j] = (float)0.0;

  1313 					break;

  1314 				case NC_DOUBLE:

  1315 					dvals[i][j] = (double)0.0;

  1316 					break;

  1317 				default:

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

  1319 					exit(3);

  1320 			}

  1321 		}

  1322 		mcsec[i] = 0;

  1323 	}

  1325 	return;

  1326 }

  1328 void free_all_samples(size_t nsamples, void **vals, int *mcsec)

  1329 {

  1330 	int i;

  1332 	for (i = 0; i < nsamples; i++)

  1333 		if (vals[i]) free(vals[i]);

  1335 	free(vals);

  1336 	free(mcsec);

  1338 	return;

  1339 }

  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)

  1342 {

  1343 	int i, ts;

  1344 	int varid;

  1345 	size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

  1346 	char name[NC_MAX_NAME+1];

  1348 	for (i = 0; i < nifnames; i++) {

  1349 #ifdef DEBUG

  1350 		printf("\tOpening %s", input_fnames[i]);

  1351 		if (ndims > 3) printf(" to retrieve level %d\n", level);

  1352 #endif /* DEBUG */

  1353 		/* Open input file */

  1354 		wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));

  1355 		/*

  1356 		 * Inquire about number of dimensions, variables, global

  1357 		 * attributes and the ID of the unlimited dimension

  1358 		 */

  1359 		wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,

  1360 			&input_ngatts, &input_unlimdimid));

  1361 		wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));

  1362 		if (strcmp(name, time_name)) {

  1363 			fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);

  1364 			exit(4);

  1365 		}

  1366 		/* Make sure we don't run off the end of allocated space */

  1367 		if ((nslice+len) > nsamples) {

  1368 			fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);

  1369 			exit(4);

  1370 		}

  1371 		/* Get seconds of day variable */

  1372 		wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));

  1373 		wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));

  1374 		/* Get variable ID for requested variable */

  1375 		wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));

  1376 		/* Retrieve time slice of desired variable */

  1377 		for (ts = 0; ts < len; ts++) {

  1378 			start[0] = ts;

  1379 			count[0] = 1;

  1380 			start[(ndims-2)] = 0;

  1381 			count[(ndims-2)] = d1;

  1382 			start[(ndims-1)] = 0;

  1383 			count[(ndims-1)] = d2;

  1384 			if (ndims == 4) {

  1385 				start[1] = level;

  1386 				count[1] = 1;

  1387 			}

  1388 			switch(nctype) {

  1389 				case NC_FLOAT:

  1390 					wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));

  1391 					break;

  1392 				case NC_DOUBLE:

  1393 					wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));

  1394 					break;

  1395 				default:

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

  1397 					exit(3);

  1398 			}

  1399 			nslice++;

  1400 		}

  1402 		/* Close input file */

  1403 		wrap_nc(nc_close(input_ncid));

  1404 	}

  1406 	return;

  1407 }

  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)

  1410 {

  1411 	int i;

  1413 	for (i = 0; i < nsamples; i++) {

  1414 		switch(nctype) {

  1415 			case NC_FLOAT:

  1416 				if (stat_type == MEAN_TYPE)

  1417 					add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);

  1418 				else

  1419 					add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);

  1420 				break;

  1421 			case NC_DOUBLE:

  1422 				if (stat_type == MEAN_TYPE)

  1423 					add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);

  1424 				else

  1425 					add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);

  1426 				break;

  1427 			default:

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

  1429 				exit(3);

  1430 		}

  1431 	}

  1433 	/* Divide sums by number of samples */

  1434 	if (stat_type == MEAN_TYPE)

  1435 		divide_means(d1, d2, means, cell_samples);

  1436 	else

  1437 		divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);

  1439 	return;

  1440 }

  1442 float *double_to_float(size_t len, double *dvar)

  1443 {

  1444 	int i;

  1445 	float *fvar;

  1447 	if (!(fvar = (float *)malloc(len * sizeof(float)))) {

  1448 		perror("double_to_float:fvar");

  1449 		exit(2);

  1450 	}

  1452 	for (i = 0; i < len; i++)

  1453 		fvar[i] = (float)dvar[i];

  1455 	return fvar;

  1456 }

  1458 void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,

  1459 	int level, int ndims, double **var_hours)

  1460 {

  1461 	/* Output dimensions should be one larger than input dimensions */

  1462 	int i, hr;

  1463 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

  1464 	float *fvar = NULL;

  1466 	if (nctype == NC_FLOAT) {

  1467 		if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {

  1468 			perror("write_var_hours:fvar");

  1469 			exit(2);

  1470 		}

  1471 	}

  1473 	for (hr = 0; hr < HOURS_PER_DAY; hr++) {

  1474 		start[0] = 0;

  1475 		count[0] = 1; /* time */

  1476 		start[1] = hr;

  1477 		count[1] = 1; /* hour */

  1478 		start[(ndims-2)] = 0;

  1479 		count[(ndims-2)] = d1;

  1480 		start[(ndims-1)] = 0;

  1481 		count[(ndims-1)] = d2;

  1482 		if (ndims == 5) {

  1483 			start[2] = level;

  1484 			count[2] = 1;

  1485 		}

  1486 		switch (nctype) {

  1487 			case NC_FLOAT:

  1488 				for (i = 0; i < (d1 * d2); i++)

  1489 					fvar[i] = (float)var_hours[hr][i];

  1490 				wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));

  1491 				break;

  1492 			case NC_DOUBLE:

  1493 				wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));

  1494 				break;

  1495 			default:

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

  1497 				exit(3);

  1498 		}

  1499 	}

  1501 	if (nctype == NC_FLOAT)

  1502 		free(fvar);

  1504 	return;

  1505 }

  1507 void compute_stats(int nifnames, char **input_fnames, size_t nsamples)

  1508 {

  1509 	int j, k, nlevels, *mcsec;

  1510 	size_t d1, d2, **cell_samples;

  1511 	double **means, **stddevs;

  1512 	struct var *in_vnode, *out_vnode;

  1513 	void **vals;

  1515 	if (input_var_head) {

  1516 		in_vnode = input_var_head;

  1517 		/* March through non-metadata variables to compute stats */

  1518 		for (j = 0; j == 0 || in_vnode != input_var_head; j++) {

  1519 			if (!in_vnode->metadata) {

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

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

  1522 					fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);

  1523 					exit(3);

  1524 				}

  1525 				/* Figure out input dimensions */

  1526 				if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {

  1527 					fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);

  1528 					exit(3);

  1529 				}

  1530 				/* Assume 2D output; find dimensions */

  1531 				d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;

  1532 				d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;

  1533 				if (in_vnode->ndims == 3)

  1534 					nlevels = 1;

  1535 				else

  1536 					nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;

  1537 				/* Allocate working space for means and stddevs */

  1538 				alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);

  1539 				init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);

  1540 				/* Allocate and initize space for entire field across all files */

  1541 				alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);

  1542 				init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);

  1543 				printf("Computing statistics for %s\n",

  1544 					in_vnode->name);

  1545 				/* Compute means and stddevs, then write them */

  1546 				for (k = 0; k < nlevels; k++) {

  1547 					/* Read the entire fields from all files */

  1548 					read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);

  1549 					/* Compute means */

  1550 					compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);

  1551 					/* Find corresponding output variable on the mean output file */

  1552 					out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);

  1553 					/* Write out the means for this variable */

  1554 					write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);

  1555 					/* Zero out cell_samples so they can be used as a flag again for computing stddevs */

  1556 					reset_cell_samples(d1, d2, cell_samples);

  1557 					/* Compute stddevs using means */

  1558 					compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);

  1559 					/* Find corresponding output variable on the stddev output file */

  1560 					out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);

  1561 					/* Write out stddevs for this variable */

  1562 					write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);

  1563 				}

  1565 				/* Free all samples */

  1566 				free_all_samples(nsamples, vals, mcsec);

  1567 				/* Free working space for means and stddevs */

  1568 				free_means_stddevs(means, stddevs, cell_samples);

  1569 			}

  1570 			in_vnode = in_vnode->next;

  1571 		}

  1572 	}

  1573 	return;

  1574 }

  1576 void usage(char *program)

  1577 {

  1578 	fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);

  1579 	fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");

  1580 	exit(4);

  1581 }

  1583 int main(int argc, char **argv)

  1584 {

  1585 	int i, ic, nifnames;

  1586 	size_t len, pos, nsamples;

  1587 	char *mfname, *sfname, **ifnames, *flist;

  1589 	ifnames = NULL;

  1590 	mfname = sfname = flist = NULL;

  1591 	nifnames = 0;

  1593 #ifdef DEBUG

  1594 	printf("Number of metadata variables (nmvars) = %d\n", nmvars);

  1595 	fflush(stdout);

  1596 #endif /* DEBUG */

  1598 	/* check command-line flags and switches */

  1599 	while ((ic = getopt(argc, argv, "m:s:")) != -1) {

  1600 		switch(ic) {

  1601 			case 'm':

  1602 				if (!(mfname = strdup(optarg))) {

  1603 					perror("mfname");

  1604 					exit(2);

  1605 				}

  1606 				break;

  1607 			case 's':

  1608 				if (!(sfname = strdup(optarg))) {

  1609 					perror("sfname");

  1610 					exit(2);

  1611 				}

  1612 				break;

  1613 		}

  1614 	}

  1616 	if (!mfname) {

  1617 		fprintf(stderr, "Output file name for writing means required.\n");

  1618 		usage(argv[0]);

  1619 	}

  1620 	if (!sfname) {

  1621 		fprintf(stderr, "Output file name for writing standard deviations required.\n");

  1622 		usage(argv[0]);

  1623 	}

  1625 	if (optind < argc) {

  1626 		if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {

  1627 			perror("ifnames");

  1628 			exit(2);

  1629 		}

  1630 		for (i = optind; i < argc; i++) {

  1631 			if (!(ifnames[nifnames++] = strdup(argv[i]))) {

  1632 				perror("ifnames[nifnames++]");

  1633 				exit(2);

  1634 			}

  1635 		}

  1636 	}

  1637 	else {

  1638 		fprintf(stderr, "At least one input file name is required.\n");

  1639 		usage(argv[0]);

  1640 	}

  1643 	/*

  1644 	 * Build list of input files to be included in the global history

  1645 	 * attribute on the two outputs files.

  1646 	 */

  1647 	for (i = len = 0; i < nifnames; i++)

  1648 		len += strlen(ifnames[i]);

  1649 	len += nifnames + 1; /* space in front of every name + null terminator */

  1650 	if (!(flist = (char *)malloc(len * sizeof(char)))) {

  1651 		perror("flist");

  1652 		exit(2);

  1653 	}

  1654 	for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)

  1655 		sprintf(&flist[pos], " %s", ifnames[i]);

  1656 #ifdef DEBUG

  1657 	printf("flist=%s\n", flist);

  1658 	fflush(stdout);

  1659 #endif /* DEBUG */

  1661 	/* Open every file to count up the number of time samples. */

  1662 	nsamples = count_samples(nifnames, ifnames);

  1663 #ifdef DEBUG

  1664 	printf("Number of samples across all files = %ld\n", (long)nsamples);

  1665 	fflush(stdout);

  1666 #endif /* DEBUG */

  1668 	/*

  1669 	 * Open the *last* input file and the two output files (for mean and

  1670 	 * standard deviation).  Define dimensions, variables, and attributes

  1671 	 * for the two output files based on the *last* input files.  The

  1672 	 * last file is used because some metadata variables will contain

  1673 	 * only values from the last time slice from the period over which

  1674 	 * the statistics are computed.

  1675 	 */

  1676 	open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,

  1677 		nsamples);

  1679 	compute_stats(nifnames, ifnames, nsamples);

  1681 	/* Close the two output files */

  1682 	wrap_nc(nc_close(mean_ncid));

  1683 	wrap_nc(nc_close(stddev_ncid));

  1685 	return 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--