c-lamp_postproc: h1_summary2.c@3c02cce30be8

     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 *mcsec_name = "mcsec";

    86 static char *history_name = "history";

    87 static char *nsamples_name = "nsamples";

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

    89 	*hour_units = "hour";

    90 static char *cell_method_name = "cell_method";

    91 /* input stuff */

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

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

    94 static struct text_att *input_gatt_head = NULL;

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

    96 static struct var *input_var_head = NULL;

    97 /* translation stuff */

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

    99 /* output stuff */

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

   101 static int mean_hour_dimid; /* special notes */

   102 static struct dim *mean_dim_head = NULL;

   103 static struct var *mean_var_head = NULL;

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

   105 static int stddev_hour_dimid; /* special notes */

   106 static struct dim *stddev_dim_head = NULL;

   107 static struct var *stddev_var_head = NULL;

   109 char is_metadata(char *name)

   110 {

   111 	int i;

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

   115 	if (i < nmvars)

   116 		return 1;

   117 	else

   118 		return 0;

   119 }

   121 #if 0

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

   123 {

   124 	int i;

   125 	struct dim *p = NULL;

   127 	if (head) {

   128 		node = head;

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

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

   131 				p = node;

   132 			node = node->next;

   133 		}

   134 		return p;

   135 	}

   136 	else

   137 		return NULL;

   138 }

   139 #endif

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

   142 {

   143 	int i;

   144 	struct var *node;

   146 	if (head) {

   147 		node = head;

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

   149 			node = node->next;

   150 		if (i && node == head)

   151 			return NULL;

   152 		else

   153 			return node;

   154 	}

   155 	else

   156 		return NULL;

   157 }

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

   160 {

   161 	struct text_att *head, *node;

   163 	head = *headp;

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

   166 		perror("gattlist_add_node");

   167 		exit(2);

   168 	}

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

   171 		perror("gattlist_add_node");

   172 		exit(2);

   173 	}

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

   175 		perror("gattlist_add_node");

   176 		exit(2);

   177 	}

   179 	if (head == NULL) {

   180 		/* first one */

   181 		*headp = node;

   182 		node->next = node;

   183 		node->prev = node;

   184 	}

   185 	else {

   186 		node->next = head;

   187 		node->prev = head->prev;

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

   189 		head->prev = node;

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

   191 		node->prev->next = node;

   192 	}

   194 	return;

   195 }

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

   198 {

   199 	struct dim *head, *node;

   201 	head = *headp;

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

   204 		perror("dimlist_add_node");

   205 		exit(2);

   206 	}

   208 	node->dimid = dimid;

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

   210 		perror("dimlist_add_node");

   211 		exit(2);

   212 	}

   213 	node->len = len;

   215 	if (head == NULL) {

   216 		/* first one */

   217 		*headp = node;

   218 		node->next = node;

   219 		node->prev = node;

   220 	}

   221 	else {

   222 		node->next = head;

   223 		node->prev = head->prev;

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

   225 		head->prev = node;

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

   227 		node->prev->next = node;

   228 	}

   230 	return node;

   231 }

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

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

   235 	float FillValue)

   236 {

   237 	int i;

   238 	struct var *head, *node;

   240 	head = *headp;

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

   243 		perror("varlist_add_node");

   244 		exit(3);

   245 	}

   247 	node->ncvarid = ncvarid;

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

   249 		perror("varlist_add_node");

   250 		exit(3);

   251 	}

   252 	node->nctype = nctype;

   253 	node->ndims = ndims;

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

   255 		perror("varlist_add_node");

   256 		exit(3);

   257 	}

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

   259 	node->natts = natts;

   260 	node->metadata = is_metadata(name);

   261 	node->FillFlag = FillFlag;

   262 	node->FillValue = FillValue;

   264 	if (head == NULL) {

   265 		/* first one */

   266 		*headp = node;

   267 		node->next = node;

   268 		node->prev = node;

   269 	}

   270 	else {

   271 		node->next = head;

   272 		node->prev = head->prev;

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

   274 		head->prev = node;

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

   276 		node->prev->next = node;

   277 	}

   279 	return;

   280 }

   282 void varlist_print(struct var *headp)

   283 {

   284 	int i, j;

   285 	struct var *node;

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

   289 	if (headp) {

   290 		node = headp;

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

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

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

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

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

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

   297 			printf("\tdimids =");

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

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

   300 			printf("\n");

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

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

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

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

   305 			node = node->next;

   306 		}

   307 	}

   308 	else {

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

   310 	}

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

   314 	return;

   315 }

   317 void wrap_nc(int status)

   318 {

   319 	if (status != NC_NOERR) {

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

   321 		exit(1);

   322 	}

   324 	return;

   325 }

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

   328 {

   329 	int i;

   330 	char name[NC_MAX_NAME+1];

   331 	size_t len;

   332 	struct dim **in_dim_idx;

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

   335 		perror("*in_dim_idxp");

   336 		exit(3);

   337 	}

   338 	in_dim_idx = *in_dim_idxp;

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

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

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

   343 	}

   345 	return;

   346 }

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

   348 {

   349 	int i;

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

   351 	nc_type xtype;

   352 	size_t len;

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

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

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

   357 		if (xtype != NC_CHAR) {

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

   359 			exit(2);

   360 		}

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

   362 			perror("get_gatts: value");

   363 			exit(3);

   364 		}

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

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

   367 		gattlist_add_node(gatt_headp, name, value);

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

   369 	}

   371 	return;

   372 }

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

   375 {

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

   377 	size_t len;

   378 	float FillValue;

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

   380 	nc_type xtype, atype;

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

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

   384 			&natts));

   385 		/* Check for _FillValue */

   386 		FillFlag = 0;

   387 		FillValue = 0.;

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

   389 			== NC_NOERR) {

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

   391 				wrap_nc(nc_get_att_float(ncid, i,

   392 					fill_att_name, &FillValue));

   393 				FillFlag = 1;

   394 			}

   395 		}

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

   397 			natts, FillFlag, FillValue);

   398 	}

   400 	return;

   401 }

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

   404 	size_t nsamples, int **in2outp, int out_ncid,

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

   406 {

   407 	/*

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

   409 	 * dimension IDs.

   410 	 */

   411 	int j, dimid, ndims, *in2out;

   412 	size_t len;

   413 	struct dim *dnode;

   415 	ndims = 0;

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

   418 		perror("put_dimensions: in2outp");

   419 		exit(3);

   420 	}

   421 	in2out = *in2outp;

   423 	if (in_dim_head) {

   424 		dnode = in_dim_head;

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

   426 			if (dnode->dimid == in_unlimdimid)

   427 				len = NC_UNLIMITED;

   428 			else

   429 				len = dnode->len;

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

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

   432 			++ndims;

   433 			in2out[dnode->dimid] = dimid;

   434 			if (dnode->dimid == in_unlimdimid)

   435 				*out_unlimdimidp = dimid;

   436 			/*

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

   438 			 * "nsamples" and "hour" dimensions.

   439 			 */

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

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

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

   443 				++ndims;

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

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

   447 				++ndims;

   448 				/* Track hour dimid for out files */

   449 				*out_hour_dimidp = dimid;

   450 			}

   452 			dnode = dnode->next;

   453 		}

   454 	}

   455 	else {

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

   457 	}

   459 	return ndims;

   460 }

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

   463 {

   464 	/*

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

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

   467 	 */

   468 	int j, hflag = 0, ngatts;

   469 	char *value;

   470 	struct text_att *anode;

   472 	ngatts = 0;

   474 	if (in_gatt_head) {

   475 		anode = in_gatt_head;

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

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

   478 				value = out_history;

   479 				++hflag;

   480 			}

   481 			else

   482 				value = anode->value;

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

   484 			++ngatts;

   485 			anode = anode->next;

   486 		}

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

   488 		if (!hflag) {

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

   490 			++ngatts;

   491 		}

   492 	}

   493 	else {

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

   495 	}

   497 	return ngatts;

   498 }

   500 int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)

   501 {

   502 	/*

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

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

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

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

   507 	 */

   508 	int i, ndims;

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

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

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

   513 			&& !metadata) {

   514 			ndims++;

   515 			out_dimids[ndims] = hour_dimid;

   516 		}

   517 	}

   519 	return ndims;

   520 }

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

   523 	int in_varid, int out_ncid, int out_varid)

   524 {

   525 	/*

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

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

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

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

   530 	 */

   532 	int i, natts;

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

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

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

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

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

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

   540 			cmflag = 1;

   541 		}

   542 		else

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

   544 	}

   545 	/*

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

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

   548 	 * on the output file.

   549 	 */

   550 	if (!cmflag && !metadata) {

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

   552 		natts++;

   553 	}

   555 	return natts;

   556 }

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

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

   560 	int stat_type)

   561 {

   562 	/*

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

   564 	 * Return number of variables defined.

   565 	 */

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

   567 	struct var *vnode;

   569 	nvars = 0;

   571 	if (in_var_head) {

   572 		vnode = in_var_head;

   573 		/*

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

   575 		 * variables on the output file.

   576 		 */

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

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

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

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

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

   582 			++nvars;

   583 			/*

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

   585 			 * variable for hour of day.

   586 			 */

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

   588 				ndims = 1;

   589 				dimids[0] = hour_dimid;

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

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

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

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

   594 				++nvars;

   595 			}

   597 			vnode = vnode->next;

   598 		}

   599 	}

   600 	else {

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

   602 	}

   604 	return nvars;

   605 }

   607 void *alloc_var(nc_type nctype, size_t len)

   608 {

   609 	void *val;

   611 	switch (nctype) {

   612 		case NC_FLOAT:

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

   614 				perror("alloc_var: val");

   615 				exit(3);

   616 			}

   617 			break;

   618 		case NC_INT:

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

   620 				perror("alloc_var: val");

   621 				exit(3);

   622 			}

   623 			break;

   624 		case NC_DOUBLE:

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

   626 				perror("alloc_var: val");

   627 				exit(3);

   628 			}

   629 			break;

   630 		case NC_CHAR:

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

   632 				perror("alloc_var: val");

   633 				exit(3);

   634 			}

   635 			break;

   636 		default:

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

   638 			exit(3);

   639 	}

   641 	return val;

   642 }

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

   645 {

   646 	int i;

   647 	size_t len = (size_t)1;

   648 	void *val;

   650 	/* Figure out the total size */

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

   653 	val = alloc_var(nctype, len);

   655 	switch (nctype) {

   656 		case NC_FLOAT:

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

   658 			break;

   659 		case NC_INT:

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

   661 			break;

   662 		case NC_DOUBLE:

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

   664 			break;

   665 		case NC_CHAR:

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

   667 			break;

   668 		default:

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

   670 			exit(3);

   671 	}

   673 	return val;

   674 }

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

   677 {

   678 	int i;

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

   680 	void *val;

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

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

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

   685 		exit(3);

   686 	}

   687 	/*

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

   689 	 * and build start/count vectors.

   690 	 */

   691 	start[0] = tslice;

   692 	count[0] = 1;

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

   694 		start[i] = 0;

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

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

   697 	}

   699 	val = alloc_var(nctype, len);

   701 	switch (nctype) {

   702 		case NC_FLOAT:

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

   704 			break;

   705 		case NC_INT:

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

   707 			break;

   708 		case NC_DOUBLE:

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

   710 			break;

   711 		case NC_CHAR:

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

   713 			break;

   714 		default:

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

   716 			exit(3);

   717 	}

   719 	return val;

   720 }

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

   723 {

   724 	switch (nctype) {

   725 		case NC_FLOAT:

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

   727 			break;

   728 		case NC_INT:

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

   730 			break;

   731 		case NC_DOUBLE:

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

   733 			break;

   734 		case NC_CHAR:

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

   736 			break;

   737 		default:

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

   739 			exit(3);

   740 	}

   742 	return;

   743 }

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

   746 {

   747 	int i;

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

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

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

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

   753 		exit(3);

   754 	}

   756 	/* Build start/count vectors */

   757 	start[0] = tslice;

   758 	count[0] = 1;

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

   760 		start[i] = 0;

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

   762 	}

   764 	switch (nctype) {

   765 		case NC_FLOAT:

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

   767 			break;

   768 		case NC_INT:

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

   770 			break;

   771 		case NC_DOUBLE:

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

   773 			break;

   774 		case NC_CHAR:

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

   776 			break;

   777 		default:

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

   779 			exit(3);

   780 	}

   782 	return;

   783 }

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

   786 	struct dim **in_dim_idx, int out_ncid, struct var *out_var_head)

   787 {

   788 	int i, j;

   789 	float hr[HOURS_PER_DAY];

   790 	struct var *in_vnode, *out_vnode;

   791 	void *val;

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

   795 	if (in_var_head) {

   796 		in_vnode = in_var_head;

   797 		/*

   798 		 * March through input variables to stuff metadata values into

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

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

   801 		 * input files).

   802 		 */

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

   804 			if (in_vnode->metadata) {

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

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

   807 					/* time is not the first dimension */

   808 #ifdef DEBUG

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

   810 						in_vnode->name);

   811 #endif /* DEBUG */

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

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

   814 				}

   815 				else {

   816 					/* time is the first dimension */

   817 #ifdef DEBUG

   818 					printf("Copying last value of \

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

   820 #endif /* DEBUG */

   821 					val = read_timeslice(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, ((input_dim_idx[in_vnode->dimids[0]])->len - 1));

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

   823 				}

   824 				free(val);

   825 				/*

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

   827 				 * the "hour" variable values.

   828 				 */

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

   830 					out_vnode = varlist_find_by_name(out_var_head, hour_name);

   831 					write_var(out_ncid, out_vnode->ncvarid,

   832 						out_vnode->nctype, hr);

   833 				}

   834 			}

   835 			else {

   836 #ifdef DEBUG

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

   838 					in_vnode->name);

   839 #endif /* DEBUG */

   840 			}

   841 			in_vnode = in_vnode->next;

   842 		}

   843 	}

   845 	return;

   846 }

   848 void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)

   849 {

   850 	char *mean_history_gatt, *stddev_history_gatt,

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

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

   854 	/*

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

   856 	 * files.

   857 	 */

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

   859 		perror("open_inout:mean_history_gatt");

   860 		exit(2);

   861 	}

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

   863 		perror("open_inout:stddev_history_gatt");

   864 		exit(2);

   865 	}

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

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

   869 	/* Open input file */

   870 	wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid));

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

   872 	 * and the ID of the unlimited dimension

   873 	 */

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

   875 		&input_unlimdimid));

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

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

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

   879 	get_gatts(input_ncid, input_ngatts, &input_gatt_head);

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

   881 	get_vars(input_ncid, input_nvars, &input_var_head);

   882 #ifdef DEBUG

   883 	varlist_print(input_var_head);

   884 #endif /* DEBUG */

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

   886 	/* Create new netCDF files */

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

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

   889 	/* Define dimensions */

   890 	mean_ndims = put_dimensions(input_dim_head, input_ndims,

   891 		input_unlimdimid, nsamples, &idid2mdid, mean_ncid,

   892 		&mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);

   893 	stddev_ndims = put_dimensions(input_dim_head, input_ndims,

   894 		input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,

   895 		&stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);

   896 	/* Define variables and their attributes */

   897 	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   898 		idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,

   899 		MEAN_TYPE);

   900 	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

   901 		idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,

   902 		STDDEV_TYPE);

   903 	/* Store global attributes */

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

   905 	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,

   906 		stddev_history_gatt);

   907 	/* End define mode */

   908 	wrap_nc(nc_enddef(mean_ncid));

   909 	wrap_nc(nc_enddef(stddev_ncid));

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

   911 	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,

   912 		mean_var_head);

   913 	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,

   914 		stddev_var_head);

   916 	wrap_nc(nc_close(input_ncid));

   918 	return;

   919 }

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

   922 {

   923 	int i;

   924 	char name[NC_MAX_NAME+1];

   925 	size_t len, cnt = 0;

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

   928 		/* Open input file */

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

   930 		/*

   931 		 * Inquire about number of dimensions, variables, global

   932 		 * attributes and the ID of the unlimited dimension

   933 		 */

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

   935 			&input_ngatts, &input_unlimdimid));

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

   937 		if (strcmp(name, time_name)) {

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

   939 			exit(4);

   940 		}

   941 #ifdef DEBUG

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

   943 #endif /* DEBUG */

   944 		wrap_nc(nc_close(input_ncid));

   945 		cnt += len;

   946 	}

   947 	return cnt;

   948 }

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

   951 {

   952 	/*

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

   954 	 * hour of day.

   955 	 */

   956 	int i;

   957 	size_t **cell_samples;

   958 	double **means, **stddevs;

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

   961 		perror("alloc_means_stddevs:*meansp");

   962 		exit(2);

   963 	}

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

   965 		perror("alloc_means_stddevs:*stddevsp");

   966 		exit(2);

   967 	}

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

   969 		perror("alloc_means_stddevs:*cell_samplesp");

   970 		exit(2);

   971 	}

   973 	means = *meansp;

   974 	stddevs = *stddevsp;

   975 	cell_samples = *cell_samplesp;

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

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

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

   980 			exit(2);

   981 		}

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

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

   984 			exit(2);

   985 		}

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

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

   988 			exit(2);

   989 		}

   990 	}

   992 	return;

   993 }

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

   996 {

   997 	int hr, i, j, pos;

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

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

  1001 #pragma _CRI concurrent

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

  1003 				pos = i*d2+j;

  1004 				means[hr][pos] = FillValue;

  1005 				stddevs[hr][pos] = FillValue;

  1006 				cell_samples[hr][pos] = 0;

  1007 			}

  1008 		}

  1009 	}

  1011 	return;

  1012 }

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

  1015 {

  1016 	int i, j, hr, pos;

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

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

  1020 #pragma _CRI concurrent

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

  1022 				pos = i*d2+j;

  1023 				cell_samples[hr][pos] = 0;

  1024 			}

  1025 		}

  1026 	}

  1028 	return;

  1029 }

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

  1032 {

  1033 	/*

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

  1035 	 * hour of day.

  1036 	 */

  1037 	int i;

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

  1040 		free(means[i]);

  1041 		free(stddevs[i]);

  1042 		free(cell_samples[i]);

  1043 	}

  1045 	free(means);

  1046 	free(stddevs);

  1047 	free(cell_samples);

  1049 	return;

  1050 }

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

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

  1054 {

  1055 	int i, j, hr, pos;

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

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

  1060 #pragma _CRI concurrent

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

  1062 			pos = i*d2+j;

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

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

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

  1066 				else

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

  1068 				++cell_samples[hr][pos];

  1069 			}

  1070 		}

  1071 	}

  1073 	return;

  1074 }

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

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

  1078 {

  1079 	int i, j, hr, pos;

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

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

  1084 #pragma _CRI concurrent

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

  1086 			pos = i*d2+j;

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

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

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

  1090 				else

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

  1092 				++cell_samples[hr][pos];

  1093 			}

  1094 		}

  1095 	}

  1097 	return;

  1098 }

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

  1102 {

  1103 	int i, j, hr, pos;

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

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

  1107 #pragma _CRI concurrent

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

  1109 				pos = i*d2+j;

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

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

  1112 				}

  1113 			}

  1114 		}

  1115 	}

  1117 	return;

  1118 }

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

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

  1122 	size_t **cell_samples)

  1123 {

  1124 	int i, j, hr, pos;

  1125 	double delta;

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

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

  1130 #pragma _CRI concurrent

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

  1132 			pos = i*d2+j;

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

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

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

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

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

  1138 				else

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

  1140 				++cell_samples[hr][pos];

  1141 			}

  1142 		}

  1143 	}

  1145 	return;

  1146 }

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

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

  1150 	size_t **cell_samples)

  1151 {

  1152 	int i, j, hr, pos;

  1153 	double delta;

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

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

  1158 #pragma _CRI concurrent

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

  1160 			pos = i*d2+j;

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

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

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

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

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

  1166 				else

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

  1168 				++cell_samples[hr][pos];

  1169 			}

  1170 		}

  1171 	}

  1173 	return;

  1174 }

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

  1178 {

  1179 	int i, j, hr, pos;

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

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

  1183 #pragma _CRI concurrent

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

  1185 				pos = i*d2+j;

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

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

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

  1189 				}

  1190 			}

  1191 		}

  1192 	}

  1194 	return;

  1195 }

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

  1198 {

  1199 	int i;

  1200 	void **vals;

  1202 	/* Allocate space for all timeslices */

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

  1204 		perror("alloc_all_samples:*valsp");

  1205 		exit(2);

  1206 	}

  1207 	vals = *valsp;

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

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

  1210 	}

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

  1212 		perror("alloc_all_samples:*mcsecp");

  1213 		exit(2);

  1214 	}

  1216 	return;

  1217 }

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

  1220 {

  1221 	int i, j;

  1222         float **fvals = NULL;

  1223 	double **dvals = NULL;

  1225 	switch(nctype) {

  1226 		case NC_FLOAT:

  1227 			fvals = (float **)vals;

  1228 			break;

  1229 		case NC_DOUBLE:

  1230 			dvals = (double **)vals;

  1231 			break;

  1232 		default:

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

  1234 			exit(3);

  1235 	}

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

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

  1239 			switch(nctype) {

  1240 				case NC_FLOAT:

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

  1242 					break;

  1243 				case NC_DOUBLE:

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

  1245 					break;

  1246 				default:

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

  1248 					exit(3);

  1249 			}

  1250 		}

  1251 		mcsec[i] = 0;

  1252 	}

  1254 	return;

  1255 }

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

  1258 {

  1259 	int i;

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

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

  1264 	free(vals);

  1265 	free(mcsec);

  1267 	return;

  1268 }

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

  1271 {

  1272 	int i, ts;

  1273 	int varid;

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

  1275 	char name[NC_MAX_NAME+1];

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

  1278 #ifdef DEBUG

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

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

  1281 #endif /* DEBUG */

  1282 		/* Open input file */

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

  1284 		/*

  1285 		 * Inquire about number of dimensions, variables, global

  1286 		 * attributes and the ID of the unlimited dimension

  1287 		 */

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

  1289 			&input_ngatts, &input_unlimdimid));

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

  1291 		if (strcmp(name, time_name)) {

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

  1293 			exit(4);

  1294 		}

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

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

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

  1298 			exit(4);

  1299 		}

  1300 		/* Get seconds of day variable */

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

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

  1303 		/* Get variable ID for requested variable */

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

  1305 		/* Retrieve time slice of desired variable */

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

  1307 			start[0] = ts;

  1308 			count[0] = 1;

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

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

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

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

  1313 			if (ndims == 4) {

  1314 				start[1] = level;

  1315 				count[1] = 1;

  1316 			}

  1317 			switch(nctype) {

  1318 				case NC_FLOAT:

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

  1320 					break;

  1321 				case NC_DOUBLE:

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

  1323 					break;

  1324 				default:

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

  1326 					exit(3);

  1327 			}

  1328 			nslice++;

  1329 		}

  1331 		/* Close input file */

  1332 		wrap_nc(nc_close(input_ncid));

  1333 	}

  1335 	return;

  1336 }

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

  1339 {

  1340 	int i;

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

  1343 		switch(nctype) {

  1344 			case NC_FLOAT:

  1345 				if (stat_type == MEAN_TYPE)

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

  1347 				else

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

  1349 				break;

  1350 			case NC_DOUBLE:

  1351 				if (stat_type == MEAN_TYPE)

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

  1353 				else

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

  1355 				break;

  1356 			default:

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

  1358 				exit(3);

  1359 		}

  1360 	}

  1362 	/* Divide sums by number of samples */

  1363 	if (stat_type == MEAN_TYPE)

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

  1365 	else

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

  1368 	return;

  1369 }

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

  1372 {

  1373 	int i;

  1374 	float *fvar;

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

  1377 		perror("double_to_float:fvar");

  1378 		exit(2);

  1379 	}

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

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

  1384 	return fvar;

  1385 }

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

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

  1389 {

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

  1391 	int i, hr;

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

  1393 	float *fvar = NULL;

  1395 	if (nctype == NC_FLOAT) {

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

  1397 			perror("write_var_hours:fvar");

  1398 			exit(2);

  1399 		}

  1400 	}

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

  1403 		start[0] = 0;

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

  1405 		start[1] = hr;

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

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

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

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

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

  1411 		if (ndims == 5) {

  1412 			start[2] = level;

  1413 			count[2] = 1;

  1414 		}

  1415 		switch (nctype) {

  1416 			case NC_FLOAT:

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

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

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

  1420 				break;

  1421 			case NC_DOUBLE:

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

  1423 				break;

  1424 			default:

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

  1426 				exit(3);

  1427 		}

  1428 	}

  1430 	if (nctype == NC_FLOAT)

  1431 		free(fvar);

  1433 	return;

  1434 }

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

  1437 {

  1438 	int j, k, nlevels, *mcsec;

  1439 	size_t d1, d2, **cell_samples;

  1440 	double **means, **stddevs;

  1441 	struct var *in_vnode, *out_vnode;

  1442 	void **vals;

  1444 	if (input_var_head) {

  1445 		in_vnode = input_var_head;

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

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

  1448 			if (!in_vnode->metadata) {

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

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

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

  1452 					exit(3);

  1453 				}

  1454 				/* Figure out input dimensions */

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

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

  1457 					exit(3);

  1458 				}

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

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

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

  1462 				if (in_vnode->ndims == 3)

  1463 					nlevels = 1;

  1464 				else

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

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

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

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

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

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

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

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

  1473 					in_vnode->name);

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

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

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

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

  1478 					/* Compute means */

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

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

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

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

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

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

  1485 					reset_cell_samples(d1, d2, cell_samples);

  1486 					/* Compute stddevs using means */

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

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

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

  1490 					/* Write out stddevs for this variable */

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

  1492 				}

  1494 				/* Free all samples */

  1495 				free_all_samples(nsamples, vals, mcsec);

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

  1497 				free_means_stddevs(means, stddevs, cell_samples);

  1498 			}

  1499 			in_vnode = in_vnode->next;

  1500 		}

  1501 	}

  1502 	return;

  1503 }

  1505 void usage(char *program)

  1506 {

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

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

  1509 	exit(4);

  1510 }

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

  1513 {

  1514 	int i, ic, nifnames;

  1515 	size_t len, pos, nsamples;

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

  1518 	ifnames = NULL;

  1519 	mfname = sfname = flist = NULL;

  1520 	nifnames = 0;

  1522 #ifdef DEBUG

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

  1524 	fflush(stdout);

  1525 #endif /* DEBUG */

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

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

  1529 		switch(ic) {

  1530 			case 'm':

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

  1532 					perror("mfname");

  1533 					exit(2);

  1534 				}

  1535 				break;

  1536 			case 's':

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

  1538 					perror("sfname");

  1539 					exit(2);

  1540 				}

  1541 				break;

  1542 		}

  1543 	}

  1545 	if (!mfname) {

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

  1547 		usage(argv[0]);

  1548 	}

  1549 	if (!sfname) {

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

  1551 		usage(argv[0]);

  1552 	}

  1554 	if (optind < argc) {

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

  1556 			perror("ifnames");

  1557 			exit(2);

  1558 		}

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

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

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

  1562 				exit(2);

  1563 			}

  1564 		}

  1565 	}

  1566 	else {

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

  1568 		usage(argv[0]);

  1569 	}

  1572 	/*

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

  1574 	 * attribute on the two outputs files.

  1575 	 */

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

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

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

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

  1580 		perror("flist");

  1581 		exit(2);

  1582 	}

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

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

  1585 #ifdef DEBUG

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

  1587 	fflush(stdout);

  1588 #endif /* DEBUG */

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

  1591 	nsamples = count_samples(nifnames, ifnames);

  1592 #ifdef DEBUG

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

  1594 	fflush(stdout);

  1595 #endif /* DEBUG */

  1597 	/*

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

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

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

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

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

  1603 	 * the statistics are computed.

  1604 	 */

  1605 	open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples);

  1607 	compute_stats(nifnames, ifnames, nsamples);

  1609 	/* Close the two output files */

  1610 	wrap_nc(nc_close(mean_ncid));

  1611 	wrap_nc(nc_close(stddev_ncid));

  1613 	return 0;

  1614 }

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