c-lamp_postproc: h1_summary2.c@978f4510987d

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

   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

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

   828 				}

   829 				free(val);

   830 				/*

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

   832 				 * the "hour" variable values.

   833 				 */

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

   835 					out_vnode = varlist_find_by_name(out_var_head, hour_name);

   836 					write_var(out_ncid, out_vnode->ncvarid,

   837 						out_vnode->nctype, hr);

   838 				}

   839 			}

   840 			else {

   841 #ifdef DEBUG

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

   843 					in_vnode->name);

   844 #endif /* DEBUG */

   845 			}

   846 			in_vnode = in_vnode->next;

   847 		}

   848 	}

   850 	return;

   851 }

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

   854 {

   855 	int time_dimid, time_bounds_varid;

   856 	size_t time_len, time_bounds_index[2];

   857 	double bnd1, bnd2;

   859 	/* Open first input file */

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

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

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

   863 	/* Get length of time dimension */

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

   865 	/* Get variable ID of time_bounds variable */

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

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

   868 	 * first timeslice of the first file */

   869 	time_bounds_index[0] = time_bounds_index[1] = 0;

   870 	/* Read the value */

   871 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,

   872 		time_bounds_index, &bnd1));

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

   874 	 * open the second one */

   875 	if (strcmp(first_fname, last_fname)) {

   876 		/* Close the first input file */

   877 		wrap_nc(nc_close(input_ncid));

   878 		/* Open the last input file */

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

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

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

   882 		/* Get length of time dimension */

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

   884 		/* Get variable ID of time_bounds variable */

   885 		wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,

   886 			&time_bounds_varid));

   887 	}

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

   889 	 * last timeslice of the last file */

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

   891 	/* Read the value */

   892 	wrap_nc(nc_get_var1_double(input_ncid, time_bounds_varid,

   893 		time_bounds_index, &bnd2));

   895 	/* Close the last input file */

   896 	wrap_nc(nc_close(input_ncid));

   898 	tbounds[0] = bnd1;

   899 	tbounds[1] = bnd2;

   901 	return;

   902 }

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

   905 {

   906 	char *mean_history_gatt, *stddev_history_gatt,

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

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

   909 	double tbounds[2];

   911 	/*

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

   913 	 * files.

   914 	 */

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

   916 		perror("open_inout:mean_history_gatt");

   917 		exit(2);

   918 	}

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

   920 		perror("open_inout:stddev_history_gatt");

   921 		exit(2);

   922 	}

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

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

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

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

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

   929 	get_time_bounds(first_fname, last_fname, tbounds);

   930 #ifdef DEBUG

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

   932 #endif /* DEBUG */

   934 	/* Open last input file */

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

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

   937 	 * and the ID of the unlimited dimension

   938 	 */

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

   940 		&input_unlimdimid));

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

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

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

   944 	get_gatts(input_ncid, input_ngatts, &input_gatt_head);

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

   946 	get_vars(input_ncid, input_nvars, &input_var_head);

   947 #ifdef DEBUG

   948 	varlist_print(input_var_head);

   949 #endif /* DEBUG */

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

   951 	/* Create new netCDF files */

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

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

   954 	/* Define dimensions */

   955 	mean_ndims = put_dimensions(input_dim_head, input_ndims,

   956 		input_unlimdimid, nsamples, &idid2mdid, mean_ncid,

   957 		&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);

   958 	stddev_ndims = put_dimensions(input_dim_head, input_ndims,

   959 		input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,

   960 		&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);

   961 	/* Define variables and their attributes */

   962 	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   963 		idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,

   964 		MEAN_TYPE);

   965 	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   966 		idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,

   967 		STDDEV_TYPE);

   968 	/* Store global attributes */

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

   970 	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,

   971 		stddev_history_gatt);

   972 	/* End define mode */

   973 	wrap_nc(nc_enddef(mean_ncid));

   974 	wrap_nc(nc_enddef(stddev_ncid));

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

   976 	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,

   977 		mean_var_head, tbounds);

   978 	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,

   979 		stddev_var_head, tbounds);

   981 	wrap_nc(nc_close(input_ncid));

   983 	return;

   984 }

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

   987 {

   988 	int i;

   989 	char name[NC_MAX_NAME+1];

   990 	size_t len, cnt = 0;

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

   993 		/* Open input file */

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

   995 		/*

   996 		 * Inquire about number of dimensions, variables, global

   997 		 * attributes and the ID of the unlimited dimension

   998 		 */

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

  1000 			&input_ngatts, &input_unlimdimid));

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

  1002 		if (strcmp(name, time_name)) {

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

  1004 			exit(4);

  1005 		}

  1006 #ifdef DEBUG

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

  1008 #endif /* DEBUG */

  1009 		wrap_nc(nc_close(input_ncid));

  1010 		cnt += len;

  1011 	}

  1012 	return cnt;

  1013 }

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

  1016 {

  1017 	/*

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

  1019 	 * hour of day.

  1020 	 */

  1021 	int i;

  1022 	size_t **cell_samples;

  1023 	double **means, **stddevs;

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

  1026 		perror("alloc_means_stddevs:*meansp");

  1027 		exit(2);

  1028 	}

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

  1030 		perror("alloc_means_stddevs:*stddevsp");

  1031 		exit(2);

  1032 	}

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

  1034 		perror("alloc_means_stddevs:*cell_samplesp");

  1035 		exit(2);

  1036 	}

  1038 	means = *meansp;

  1039 	stddevs = *stddevsp;

  1040 	cell_samples = *cell_samplesp;

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

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

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

  1045 			exit(2);

  1046 		}

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

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

  1049 			exit(2);

  1050 		}

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

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

  1053 			exit(2);

  1054 		}

  1055 	}

  1057 	return;

  1058 }

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

  1061 {

  1062 	int hr, i, j, pos;

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

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

  1066 #pragma _CRI concurrent

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

  1068 				pos = i*d2+j;

  1069 				means[hr][pos] = FillValue;

  1070 				stddevs[hr][pos] = FillValue;

  1071 				cell_samples[hr][pos] = 0;

  1072 			}

  1073 		}

  1074 	}

  1076 	return;

  1077 }

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

  1080 {

  1081 	int i, j, hr, pos;

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

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

  1085 #pragma _CRI concurrent

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

  1087 				pos = i*d2+j;

  1088 				cell_samples[hr][pos] = 0;

  1089 			}

  1090 		}

  1091 	}

  1093 	return;

  1094 }

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

  1097 {

  1098 	/*

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

  1100 	 * hour of day.

  1101 	 */

  1102 	int i;

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

  1105 		free(means[i]);

  1106 		free(stddevs[i]);

  1107 		free(cell_samples[i]);

  1108 	}

  1110 	free(means);

  1111 	free(stddevs);

  1112 	free(cell_samples);

  1114 	return;

  1115 }

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

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

  1119 {

  1120 	int i, j, hr, pos;

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

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

  1125 #pragma _CRI concurrent

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

  1127 			pos = i*d2+j;

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

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

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

  1131 				else

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

  1133 				++cell_samples[hr][pos];

  1134 			}

  1135 		}

  1136 	}

  1138 	return;

  1139 }

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

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

  1143 {

  1144 	int i, j, hr, pos;

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

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

  1149 #pragma _CRI concurrent

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

  1151 			pos = i*d2+j;

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

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

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

  1155 				else

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

  1157 				++cell_samples[hr][pos];

  1158 			}

  1159 		}

  1160 	}

  1162 	return;

  1163 }

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

  1167 {

  1168 	int i, j, hr, pos;

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

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

  1172 #pragma _CRI concurrent

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

  1174 				pos = i*d2+j;

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

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

  1177 				}

  1178 			}

  1179 		}

  1180 	}

  1182 	return;

  1183 }

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

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

  1187 	size_t **cell_samples)

  1188 {

  1189 	int i, j, hr, pos;

  1190 	double delta;

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

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

  1195 #pragma _CRI concurrent

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

  1197 			pos = i*d2+j;

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

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

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

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

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

  1203 				else

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

  1205 				++cell_samples[hr][pos];

  1206 			}

  1207 		}

  1208 	}

  1210 	return;

  1211 }

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

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

  1215 	size_t **cell_samples)

  1216 {

  1217 	int i, j, hr, pos;

  1218 	double delta;

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

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

  1223 #pragma _CRI concurrent

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

  1225 			pos = i*d2+j;

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

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

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

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

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

  1231 				else

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

  1233 				++cell_samples[hr][pos];

  1234 			}

  1235 		}

  1236 	}

  1238 	return;

  1239 }

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

  1243 {

  1244 	int i, j, hr, pos;

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

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

  1248 #pragma _CRI concurrent

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

  1250 				pos = i*d2+j;

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

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

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

  1254 				}

  1255 			}

  1256 		}

  1257 	}

  1259 	return;

  1260 }

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

  1263 {

  1264 	int i;

  1265 	void **vals;

  1267 	/* Allocate space for all timeslices */

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

  1269 		perror("alloc_all_samples:*valsp");

  1270 		exit(2);

  1271 	}

  1272 	vals = *valsp;

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

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

  1275 	}

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

  1277 		perror("alloc_all_samples:*mcsecp");

  1278 		exit(2);

  1279 	}

  1281 	return;

  1282 }

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

  1285 {

  1286 	int i, j;

  1287         float **fvals = NULL;

  1288 	double **dvals = NULL;

  1290 	switch(nctype) {

  1291 		case NC_FLOAT:

  1292 			fvals = (float **)vals;

  1293 			break;

  1294 		case NC_DOUBLE:

  1295 			dvals = (double **)vals;

  1296 			break;

  1297 		default:

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

  1299 			exit(3);

  1300 	}

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

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

  1304 			switch(nctype) {

  1305 				case NC_FLOAT:

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

  1307 					break;

  1308 				case NC_DOUBLE:

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

  1310 					break;

  1311 				default:

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

  1313 					exit(3);

  1314 			}

  1315 		}

  1316 		mcsec[i] = 0;

  1317 	}

  1319 	return;

  1320 }

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

  1323 {

  1324 	int i;

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

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

  1329 	free(vals);

  1330 	free(mcsec);

  1332 	return;

  1333 }

  1335 void read_all_samples(int nifnames, char **input_fnames, size_t d1, size_t d2, size_t nsamples, char *var_name, nc_type nctype, int level, int ndims, void **vals, int *mcsec)

  1336 {

  1337 	int i, ts;

  1338 	int varid;

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

  1340 	char name[NC_MAX_NAME+1];

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

  1343 #ifdef DEBUG

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

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

  1346 #endif /* DEBUG */

  1347 		/* Open input file */

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

  1349 		/*

  1350 		 * Inquire about number of dimensions, variables, global

  1351 		 * attributes and the ID of the unlimited dimension

  1352 		 */

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

  1354 			&input_ngatts, &input_unlimdimid));

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

  1356 		if (strcmp(name, time_name)) {

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

  1358 			exit(4);

  1359 		}

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

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

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

  1363 			exit(4);

  1364 		}

  1365 		/* Get seconds of day variable */

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

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

  1368 		/* Get variable ID for requested variable */

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

  1370 		/* Retrieve time slice of desired variable */

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

  1372 			start[0] = ts;

  1373 			count[0] = 1;

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

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

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

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

  1378 			if (ndims == 4) {

  1379 				start[1] = level;

  1380 				count[1] = 1;

  1381 			}

  1382 			switch(nctype) {

  1383 				case NC_FLOAT:

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

  1385 					break;

  1386 				case NC_DOUBLE:

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

  1388 					break;

  1389 				default:

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

  1391 					exit(3);

  1392 			}

  1393 			nslice++;

  1394 		}

  1396 		/* Close input file */

  1397 		wrap_nc(nc_close(input_ncid));

  1398 	}

  1400 	return;

  1401 }

  1403 void compute(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples)

  1404 {

  1405 	int i;

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

  1408 		switch(nctype) {

  1409 			case NC_FLOAT:

  1410 				if (stat_type == MEAN_TYPE)

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

  1412 				else

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

  1414 				break;

  1415 			case NC_DOUBLE:

  1416 				if (stat_type == MEAN_TYPE)

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

  1418 				else

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

  1420 				break;

  1421 			default:

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

  1423 				exit(3);

  1424 		}

  1425 	}

  1427 	/* Divide sums by number of samples */

  1428 	if (stat_type == MEAN_TYPE)

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

  1430 	else

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

  1433 	return;

  1434 }

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

  1437 {

  1438 	int i;

  1439 	float *fvar;

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

  1442 		perror("double_to_float:fvar");

  1443 		exit(2);

  1444 	}

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

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

  1449 	return fvar;

  1450 }

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

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

  1454 {

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

  1456 	int i, hr;

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

  1458 	float *fvar = NULL;

  1460 	if (nctype == NC_FLOAT) {

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

  1462 			perror("write_var_hours:fvar");

  1463 			exit(2);

  1464 		}

  1465 	}

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

  1468 		start[0] = 0;

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

  1470 		start[1] = hr;

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

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

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

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

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

  1476 		if (ndims == 5) {

  1477 			start[2] = level;

  1478 			count[2] = 1;

  1479 		}

  1480 		switch (nctype) {

  1481 			case NC_FLOAT:

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

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

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

  1485 				break;

  1486 			case NC_DOUBLE:

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

  1488 				break;

  1489 			default:

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

  1491 				exit(3);

  1492 		}

  1493 	}

  1495 	if (nctype == NC_FLOAT)

  1496 		free(fvar);

  1498 	return;

  1499 }

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

  1502 {

  1503 	int j, k, nlevels, *mcsec;

  1504 	size_t d1, d2, **cell_samples;

  1505 	double **means, **stddevs;

  1506 	struct var *in_vnode, *out_vnode;

  1507 	void **vals;

  1509 	if (input_var_head) {

  1510 		in_vnode = input_var_head;

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

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

  1513 			if (!in_vnode->metadata) {

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

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

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

  1517 					exit(3);

  1518 				}

  1519 				/* Figure out input dimensions */

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

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

  1522 					exit(3);

  1523 				}

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

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

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

  1527 				if (in_vnode->ndims == 3)

  1528 					nlevels = 1;

  1529 				else

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

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

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

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

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

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

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

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

  1538 					in_vnode->name);

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

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

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

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

  1543 					/* Compute means */

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

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

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

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

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

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

  1550 					reset_cell_samples(d1, d2, cell_samples);

  1551 					/* Compute stddevs using means */

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

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

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

  1555 					/* Write out stddevs for this variable */

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

  1557 				}

  1559 				/* Free all samples */

  1560 				free_all_samples(nsamples, vals, mcsec);

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

  1562 				free_means_stddevs(means, stddevs, cell_samples);

  1563 			}

  1564 			in_vnode = in_vnode->next;

  1565 		}

  1566 	}

  1567 	return;

  1568 }

  1570 void usage(char *program)

  1571 {

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

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

  1574 	exit(4);

  1575 }

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

  1578 {

  1579 	int i, ic, nifnames;

  1580 	size_t len, pos, nsamples;

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

  1583 	ifnames = NULL;

  1584 	mfname = sfname = flist = NULL;

  1585 	nifnames = 0;

  1587 #ifdef DEBUG

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

  1589 	fflush(stdout);

  1590 #endif /* DEBUG */

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

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

  1594 		switch(ic) {

  1595 			case 'm':

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

  1597 					perror("mfname");

  1598 					exit(2);

  1599 				}

  1600 				break;

  1601 			case 's':

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

  1603 					perror("sfname");

  1604 					exit(2);

  1605 				}

  1606 				break;

  1607 		}

  1608 	}

  1610 	if (!mfname) {

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

  1612 		usage(argv[0]);

  1613 	}

  1614 	if (!sfname) {

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

  1616 		usage(argv[0]);

  1617 	}

  1619 	if (optind < argc) {

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

  1621 			perror("ifnames");

  1622 			exit(2);

  1623 		}

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

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

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

  1627 				exit(2);

  1628 			}

  1629 		}

  1630 	}

  1631 	else {

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

  1633 		usage(argv[0]);

  1634 	}

  1637 	/*

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

  1639 	 * attribute on the two outputs files.

  1640 	 */

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

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

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

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

  1645 		perror("flist");

  1646 		exit(2);

  1647 	}

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

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

  1650 #ifdef DEBUG

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

  1652 	fflush(stdout);

  1653 #endif /* DEBUG */

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

  1656 	nsamples = count_samples(nifnames, ifnames);

  1657 #ifdef DEBUG

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

  1659 	fflush(stdout);

  1660 #endif /* DEBUG */

  1662 	/*

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

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

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

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

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

  1668 	 * the statistics are computed.

  1669 	 */

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

  1671 		nsamples);

  1673 	compute_stats(nifnames, ifnames, nsamples);

  1675 	/* Close the two output files */

  1676 	wrap_nc(nc_close(mean_ncid));

  1677 	wrap_nc(nc_close(stddev_ncid));

  1679 	return 0;

  1680 }

author	Forrest Hoffman <forrest@climatemodeling.org>
	Mon, 01 Oct 2007 15:49:25 -0400
changeset 2	978f4510987d
parent 0	3c02cce30be8
child 3	d3122367777b
permissions	-rw-r--r--