c-lamp_postproc: h1_summary

     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

    16 #define NUM_TYPES	2

    18 #define HOURS_PER_DAY	24

    19 #define MIN_PER_HOUR	60

    20 #define SEC_PER_MIN	60

    21 #define SEC_PER_HOUR	(MIN_PER_HOUR * SEC_PER_MIN)

    23 static char *metadata_vars[] = {

    24 	"lon",

    25 	"lat",

    26 	"lonrof",

    27 	"latrof",

    28 	"time",

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

    30 	"levsoi",

    31 	"levlak",

    32 	"edgen",

    33 	"edgee",

    34 	"edges",

    35 	"edgew",

    36 	"longxy",

    37 	"latixy",

    38 	"area",

    39 	"areaupsc",

    40 	"topo",

    41 	"topodnsc",

    42 	"landfrac",

    43 	"landmask",

    44 	"pftmask",

    45 	"indxupsc",

    46 	"mcdate",

    47 	"mcsec",

    48 	"mdcur",

    49 	"mscur",

    50 	"nstep",

    51 	"time_bounds",

    52 	"date_written",

    53 	"time_written"

    54 };

    56 struct text_att {

    57 	char *name;

    58 	char *value;

    59 	struct text_att *next;

    60 	struct text_att *prev;

    61 };

    63 struct dim {

    64 	int dimid;

    65 	char *name;

    66 	size_t len;

    67 	struct dim *next;

    68 	struct dim *prev;

    69 };

    71 struct var {

    72 	int ncvarid;

    73 	char *name;

    74 	nc_type nctype;

    75 	int ndims;

    76 	int *dimids;

    77 	int natts;

    78 	char metadata;

    79 	char FillFlag;

    80 	float FillValue;

    81 	struct var *next;

    82 	struct var *prev;

    83 };

    85 static char *time_name = "time";

    86 static char *time_bounds_attname = "bounds";

    87 static char *time_bounds_name = "time_bounds";

    88 static char *climatology_bounds_attname = "climatology";

    89 static char *climatology_bounds_name = "climatology_bounds";

    90 static char *mcdate_name = "mcdate";

    91 static char *mcsec_name = "mcsec";

    92 static char *history_name = "history";

    93 static char *nsamples_name = "nsamples";

    94 static char *cell_method_name = "cell_method";

    95 static char *cell_methods_prefix[] = { "time: mean within days time: mean over days", "time: mean within days times: standard_deviation over days" };

    96 /* input stuff */

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

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

    99 static struct text_att *input_gatt_head = NULL;

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

   101 static struct var *input_var_head = NULL;

   102 /* translation stuff */

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

   104 /* output stuff */

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

   106 static struct dim *mean_dim_head = NULL;

   107 static struct var *mean_var_head = NULL;

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

   109 static struct dim *stddev_dim_head = NULL;

   110 static struct var *stddev_var_head = NULL;

   111 /* output values */

   112 static double *times;

   113 static double *tbounds;

   114 static int *mcdate;

   115 static int *mcsec;

   117 char is_metadata(char *name)

   118 {

   119 	int i;

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

   123 	if (i < nmvars)

   124 		return 1;

   125 	else

   126 		return 0;

   127 }

   129 #if 0

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

   131 {

   132 	int i;

   133 	struct dim *p = NULL;

   135 	if (head) {

   136 		node = head;

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

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

   139 				p = node;

   140 			node = node->next;

   141 		}

   142 		return p;

   143 	}

   144 	else

   145 		return NULL;

   146 }

   147 #endif

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

   150 {

   151 	int i;

   152 	struct var *node;

   154 	if (head) {

   155 		node = head;

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

   157 			node = node->next;

   158 		if (i && node == head)

   159 			return NULL;

   160 		else

   161 			return node;

   162 	}

   163 	else

   164 		return NULL;

   165 }

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

   168 {

   169 	struct text_att *head, *node;

   171 	head = *headp;

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

   174 		perror("gattlist_add_node");

   175 		exit(2);

   176 	}

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

   179 		perror("gattlist_add_node");

   180 		exit(2);

   181 	}

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

   183 		perror("gattlist_add_node");

   184 		exit(2);

   185 	}

   187 	if (head == NULL) {

   188 		/* first one */

   189 		*headp = node;

   190 		node->next = node;

   191 		node->prev = node;

   192 	}

   193 	else {

   194 		node->next = head;

   195 		node->prev = head->prev;

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

   197 		head->prev = node;

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

   199 		node->prev->next = node;

   200 	}

   202 	return;

   203 }

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

   206 {

   207 	struct dim *head, *node;

   209 	head = *headp;

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

   212 		perror("dimlist_add_node");

   213 		exit(2);

   214 	}

   216 	node->dimid = dimid;

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

   218 		perror("dimlist_add_node");

   219 		exit(2);

   220 	}

   221 	node->len = len;

   223 	if (head == NULL) {

   224 		/* first one */

   225 		*headp = node;

   226 		node->next = node;

   227 		node->prev = node;

   228 	}

   229 	else {

   230 		node->next = head;

   231 		node->prev = head->prev;

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

   233 		head->prev = node;

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

   235 		node->prev->next = node;

   236 	}

   238 	return node;

   239 }

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

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

   243 	float FillValue)

   244 {

   245 	int i;

   246 	struct var *head, *node;

   248 	head = *headp;

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

   251 		perror("varlist_add_node");

   252 		exit(3);

   253 	}

   255 	node->ncvarid = ncvarid;

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

   257 		perror("varlist_add_node");

   258 		exit(3);

   259 	}

   260 	node->nctype = nctype;

   261 	node->ndims = ndims;

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

   263 		perror("varlist_add_node");

   264 		exit(3);

   265 	}

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

   267 	node->natts = natts;

   268 	node->metadata = is_metadata(name);

   269 	node->FillFlag = FillFlag;

   270 	node->FillValue = FillValue;

   272 	if (head == NULL) {

   273 		/* first one */

   274 		*headp = node;

   275 		node->next = node;

   276 		node->prev = node;

   277 	}

   278 	else {

   279 		node->next = head;

   280 		node->prev = head->prev;

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

   282 		head->prev = node;

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

   284 		node->prev->next = node;

   285 	}

   287 	return;

   288 }

   290 void varlist_print(struct var *headp)

   291 {

   292 	int i, j;

   293 	struct var *node;

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

   297 	if (headp) {

   298 		node = headp;

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

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

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

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

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

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

   305 			printf("\tdimids =");

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

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

   308 			printf("\n");

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

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

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

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

   313 			node = node->next;

   314 		}

   315 	}

   316 	else {

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

   318 	}

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

   322 	return;

   323 }

   325 void wrap_nc(int status)

   326 {

   327 	if (status != NC_NOERR) {

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

   329 		exit(1);

   330 	}

   332 	return;

   333 }

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

   336 {

   337 	int i;

   338 	char name[NC_MAX_NAME+1];

   339 	size_t len;

   340 	struct dim **in_dim_idx;

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

   343 		perror("*in_dim_idxp");

   344 		exit(3);

   345 	}

   346 	in_dim_idx = *in_dim_idxp;

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

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

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

   351 	}

   353 	return;

   354 }

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

   356 {

   357 	int i;

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

   359 	nc_type xtype;

   360 	size_t len;

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

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

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

   365 		if (xtype != NC_CHAR) {

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

   367 			exit(2);

   368 		}

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

   370 			perror("get_gatts: value");

   371 			exit(3);

   372 		}

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

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

   375 		gattlist_add_node(gatt_headp, name, value);

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

   377 	}

   379 	return;

   380 }

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

   383 {

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

   385 	size_t len;

   386 	float FillValue;

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

   388 	nc_type xtype, atype;

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

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

   392 			&natts));

   393 		/* Check for _FillValue */

   394 		FillFlag = 0;

   395 		FillValue = 0.;

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

   397 			== NC_NOERR) {

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

   399 				wrap_nc(nc_get_att_float(ncid, i,

   400 					fill_att_name, &FillValue));

   401 				FillFlag = 1;

   402 			}

   403 		}

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

   405 			natts, FillFlag, FillValue);

   406 	}

   408 	return;

   409 }

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

   412 	size_t nsamples, int **in2outp, int out_ncid,

   413 	struct dim **out_dim_headp, int *out_unlimdimidp)

   414 {

   415 	/*

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

   417 	 * dimension IDs.

   418 	 */

   419 	int j, dimid, ndims, *in2out;

   420 	size_t len;

   421 	struct dim *dnode;

   423 	ndims = 0;

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

   426 		perror("put_dimensions:in2outp");

   427 		exit(3);

   428 	}

   429 	in2out = *in2outp;

   431 	if (in_dim_head) {

   432 		dnode = in_dim_head;

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

   434 			if (dnode->dimid == in_unlimdimid)

   435 				len = NC_UNLIMITED;

   436 			else

   437 				len = dnode->len;

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

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

   440 			++ndims;

   441 			in2out[dnode->dimid] = dimid;

   442 			if (dnode->dimid == in_unlimdimid)

   443 				*out_unlimdimidp = dimid;

   444 			/* Just after the "time" dimension, add the new

   445 			 * "nsamples" dimension. */

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

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

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

   449 				++ndims;

   450 			}

   451 			dnode = dnode->next;

   452 		}

   453 	}

   454 	else {

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

   456 	}

   458 	return ndims;

   459 }

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

   462 {

   463 	/*

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

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

   466 	 */

   467 	int j, hflag = 0, ngatts;

   468 	char *value;

   469 	struct text_att *anode;

   471 	ngatts = 0;

   473 	if (in_gatt_head) {

   474 		anode = in_gatt_head;

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

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

   477 				value = out_history;

   478 				++hflag;

   479 			}

   480 			else

   481 				value = anode->value;

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

   483 			++ngatts;

   484 			anode = anode->next;

   485 		}

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

   487 		if (!hflag) {

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

   489 			++ngatts;

   490 		}

   491 	}

   492 	else {

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

   494 	}

   496 	return ngatts;

   497 }

   499 int translate_dimids(int in_ndims, int *in_dimids, int *in2out, int *out_dimids)

   500 {

   501 	/*

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

   503 	 */

   504 	int i, ndims;

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

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

   509 	return ndims;

   510 }

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

   513 	int in_varid, int out_ncid, int out_varid, char **cell_methods)

   514 {

   515 	/*

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

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

   518 	 * the cell_method attribute is set correctly, depending on the output

   519 	 * file type.

   520 	 */

   522 	int i, natts;

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

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

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

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

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

   529 			cmflag = 1;

   530 		}

   531 		else if (!strcmp(name, time_bounds_attname))

   532 			/* Change time_bounds to climatology_bounds */

   533 			wrap_nc(nc_put_att_text(out_ncid, out_varid, climatology_bounds_attname, strlen(climatology_bounds_name), climatology_bounds_name));

   534 		else

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

   536 	}

   537 	/*

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

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

   540 	 * on the output file.

   541 	 */

   542 	if (!cmflag && !metadata) {

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

   544 		natts++;

   545 	}

   547 	return natts;

   548 }

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

   551 	int *in2out, int out_ncid, struct var **out_var_headp,

   552 	char **cell_methods, int stat_type)

   553 {

   554 	/*

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

   556 	 * Return number of variables defined.

   557 	 */

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

   559 	struct var *vnode;

   561 	nvars = 0;

   563 	if (in_var_head) {

   564 		vnode = in_var_head;

   565 		/*

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

   567 		 * variables on the output file.

   568 		 */

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

   570 			ndims = translate_dimids(vnode->ndims, vnode->dimids, in2out, dimids);

   571 			/* Create all normal data variables, but only the four

   572 			 * time-based metadata variables: time,

   573 			 * climatology_bounds, mcdate, and mcsec. */

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

   575 				/* Magically promote "time" variable to

   576 				 * double precision on output file */

   577 				wrap_nc(nc_def_var(out_ncid, vnode->name, NC_DOUBLE, ndims, dimids, &varid));

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

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

   580 				++nvars;

   581 				/* Create "climatology_bounds" right after

   582 				 * time and instead of "time_bounds", with no

   583 				 * attributes */

   584 				wrap_nc(nc_def_var(out_ncid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, &varid));

   585 				natts = 0;

   586 				varlist_add_node(out_var_headp, varid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);

   587 				++nvars;

   588 			}

   589 			else if (!vnode->metadata || !strcmp(vnode->name, mcdate_name) || !strcmp(vnode->name, mcsec_name) || (vnode->metadata && strcmp((in_dim_idx[vnode->dimids[0]])->name, time_name))) {

   590 				/* If it is not a metadata variable OR it is

   591 				 * mcdate OR it is mcsec OR (it is a metadata

   592 				 * variable that does not have time as its

   593 				 * first dimension), then create it */

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

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

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

   597 				++nvars;

   598 			}

   599 			vnode = vnode->next;

   600 		}

   601 	}

   602 	else {

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

   604 	}

   606 	return nvars;

   607 }

   609 void *alloc_var(nc_type nctype, size_t len)

   610 {

   611 	void *val;

   613 	switch (nctype) {

   614 		case NC_FLOAT:

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

   616 				perror("alloc_var: val");

   617 				exit(3);

   618 			}

   619 			break;

   620 		case NC_INT:

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

   622 				perror("alloc_var: val");

   623 				exit(3);

   624 			}

   625 			break;

   626 		case NC_DOUBLE:

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

   628 				perror("alloc_var: val");

   629 				exit(3);

   630 			}

   631 			break;

   632 		case NC_CHAR:

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

   634 				perror("alloc_var: val");

   635 				exit(3);

   636 			}

   637 			break;

   638 		default:

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

   640 			exit(3);

   641 	}

   643 	return val;

   644 }

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

   647 {

   648 	int i;

   649 	size_t len = (size_t)1;

   650 	void *val;

   652 	/* Figure out the total size */

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

   655 	val = alloc_var(nctype, len);

   657 	switch (nctype) {

   658 		case NC_FLOAT:

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

   660 			break;

   661 		case NC_INT:

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

   663 			break;

   664 		case NC_DOUBLE:

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

   666 			break;

   667 		case NC_CHAR:

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

   669 			break;

   670 		default:

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

   672 			exit(3);

   673 	}

   675 	return val;

   676 }

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

   679 {

   680 	int i;

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

   682 	void *val;

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

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

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

   687 		exit(3);

   688 	}

   689 	/*

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

   691 	 * and build start/count vectors.

   692 	 */

   693 	start[0] = tslice;

   694 	count[0] = 1;

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

   696 		start[i] = 0;

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

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

   699 	}

   701 	val = alloc_var(nctype, len);

   703 	switch (nctype) {

   704 		case NC_FLOAT:

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

   706 			break;

   707 		case NC_INT:

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

   709 			break;

   710 		case NC_DOUBLE:

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

   712 			break;

   713 		case NC_CHAR:

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

   715 			break;

   716 		default:

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

   718 			exit(3);

   719 	}

   721 	return val;

   722 }

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

   725 {

   726 	switch (nctype) {

   727 		case NC_FLOAT:

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

   729 			break;

   730 		case NC_INT:

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

   732 			break;

   733 		case NC_DOUBLE:

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

   735 			break;

   736 		case NC_CHAR:

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

   738 			break;

   739 		default:

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

   741 			exit(3);

   742 	}

   744 	return;

   745 }

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

   748 {

   749 	int i;

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

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

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

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

   755 		exit(3);

   756 	}

   758 	/* Build start/count vectors */

   759 	start[0] = tslice;

   760 	count[0] = 1;

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

   762 		start[i] = 0;

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

   764 	}

   766 	switch (nctype) {

   767 		case NC_FLOAT:

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

   769 			break;

   770 		case NC_INT:

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

   772 			break;

   773 		case NC_DOUBLE:

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

   775 			break;

   776 		case NC_CHAR:

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

   778 			break;

   779 		default:

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

   781 			exit(3);

   782 	}

   784 	return;

   785 }

   787 void write_time_var(int ncid, int varid, nc_type nctype, void *time_var)

   788 {

   789 	size_t start[1], count[1];

   791 	start[0] = 0;

   792 	count[0] = 24;

   794 	switch (nctype) {

   795 		case NC_FLOAT:

   796 			wrap_nc(nc_put_vara_float(ncid, varid, start, count, time_var));

   797 			break;

   798 		case NC_INT:

   799 			wrap_nc(nc_put_vara_int(ncid, varid, start, count, time_var));

   800 			break;

   801 		case NC_DOUBLE:

   802 			wrap_nc(nc_put_vara_double(ncid, varid, start, count, time_var));

   803 			break;

   804 		case NC_CHAR:

   805 			wrap_nc(nc_put_vara_text(ncid, varid, start, count, time_var));

   806 			break;

   807 		default:

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

   809 			exit(3);

   810 	}

   812 	return;

   813 }

   815 void write_tbounds(int ncid, int varid, double *tbounds)

   816 {

   817 	size_t start[2], count[2];

   819 	start[0] = start[1] = 0;

   820 	count[0] = 24; count[1] = 2;

   822 	wrap_nc(nc_put_vara_double(ncid, varid, start, count, tbounds));

   824 	return;

   825 }

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

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

   829 	double *times, double *tbounds, int *mcdate, int *mcsec)

   830 {

   831 	int i;

   832 	struct var *in_vnode, *out_vnode;

   833 	void *val;

   835 	if (!in_var_head) return;

   837 	in_vnode = in_var_head;

   838 	/*

   839 	 * March through input variables to stuff metadata values into

   840 	 * the new output files. NOTE: The only time-based metadata

   841 	 * variables, explicitly handled, are: time, climatology_bounds,

   842 	 * mcdate, and mcsec.  These are passed in.

   843 	 */

   844 	for (i = 0; i == 0 || in_vnode != in_var_head; i++) {

   845 		if (in_vnode->metadata) {

   846 			/* Find the corresponding output variable */

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

   848 				out_vnode = varlist_find_by_name(out_var_head, climatology_bounds_name);

   849 			else

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

   851 			if (out_vnode) {

   852 				/* Read values from input and write to output */

   853 				if (strcmp((in_dim_idx[in_vnode->dimids[0]])->name, time_name)) {

   854 					/* "time" is not the first dimension,

   855 					 * so just copy the variable values */

   856 #ifdef DEBUG

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

   858 						in_vnode->name);

   859 #endif /* DEBUG */

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

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

   862 					free(val);

   863 				}

   864 				else {

   865 					/* "time" is the first dimension, so

   866 					 * explicitly handle the "time",

   867 					 * "climatology_bounds", "mcdate", and

   868 					 * "mcsec" while dropping all others */

   869 					if (!strcmp(out_vnode->name, time_name)) {

   870 						/* Force time stamps to be

   871 						 * those computed previously,

   872 						 * and remember these are now

   873 						 * double precision. */

   874 #ifdef DEBUG

   875 						printf("Writing metadata variable %s\n", out_vnode->name);

   876 #endif /* DEBUG */

   877 						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, times);

   878 					}

   879 					if (!strcmp(out_vnode->name, climatology_bounds_name)) {

   880 						/* time_bounds are now

   881 						 * climatology_bounds and have

   882 						 * pre-computed values */

   883 #ifdef DEBUG

   884 						printf("Writing metadata variable %s\n", out_vnode->name);

   885 #endif /* DEBUG */

   886 						write_tbounds(out_ncid, out_vnode->ncvarid, tbounds);

   887 					}

   888 					else if (!strcmp(out_vnode->name, mcdate_name)) {

   889 						/* mcdate is pre-computed */

   890 #ifdef DEBUG

   891 						printf("Writing metadata variable %s\n", out_vnode->name);

   892 #endif /* DEBUG */

   893 						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcdate);

   894 					}

   895 					else if (!strcmp(out_vnode->name, mcsec_name)) {

   896 						/* mcsec is pre-computed */

   897 #ifdef DEBUG

   898 						printf("Writing metadata variable %s\n", out_vnode->name);

   899 #endif /* DEBUG */

   900 						write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcsec);

   901 					}

   902 				}

   903 			}

   904 			else {

   905 #ifdef DEBUG

   906 				printf("Dropping metadata variable %s\n", in_vnode->name);

   907 #endif /* DEBUG */

   908 			}

   909 		}

   910 		else {

   911 #ifdef DEBUG

   912 			printf("Skipping non-metadata variable %s\n",

   913 				in_vnode->name);

   914 #endif /* DEBUG */

   915 		}

   916 #ifdef DEBUG

   917 		printf ("Done\n");

   918 #endif /* DEBUG */

   919 		in_vnode = in_vnode->next;

   920 	}

   922 	return;

   923 }

   925 void get_time_bounds(char *first_fname, char *last_fname, double *times, double *tbounds, int *mcdate, int *mcsec)

   926 {

   927 	int i, time_dimid, time_bounds_varid, mcdate_varid, mcsec_varid;

   928 	size_t time_len, time_bounds_start[2], time_bounds_count[2],

   929 		mc_start[1], mc_count[1];

   930 	double bnd1[24], bnd2[24], day1, day2;

   932 	/* Open first input file */

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

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

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

   936 	/* Get length of time dimension */

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

   938 	/* Get variable ID of time_bounds variable */

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

   940 	/* Set start/count for reading the left values of time_bounds from the

   941 	 * first 24 timeslices of the first file */

   942 	time_bounds_start[0] = time_bounds_start[1] = 0;

   943 	time_bounds_count[0] = 24; time_bounds_count[1] = 1;

   944 	/* Read the values */

   945 	wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,

   946 		time_bounds_start, time_bounds_count, bnd1));

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

   948 	 * open the second one */

   949 	if (strcmp(first_fname, last_fname)) {

   950 		/* Close the first input file */

   951 		wrap_nc(nc_close(input_ncid));

   952 		/* Open the last input file */

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

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

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

   956 		/* Get length of time dimension */

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

   958 		/* Get variable ID of time_bounds variable */

   959 		wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,

   960 			&time_bounds_varid));

   961 	}

   962 	/* Set start/count for reading the right values of time_bounds from the

   963 	 * last 24 timeslices of the last file */

   964 	time_bounds_start[0] = time_len - 24; time_bounds_start[1] = 1;

   965 	time_bounds_count[0] = 24; time_bounds_count[1] = 1;

   966 	/* Read the value */

   967 	wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,

   968 		time_bounds_start, time_bounds_count, bnd2));

   969 	/* Read the last 24 mcdate and mcsec values */

   970 	wrap_nc(nc_inq_varid(input_ncid, mcdate_name, &mcdate_varid));

   971 	wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &mcsec_varid));

   972 	mc_start[0] = time_len - 24;

   973 	mc_count[0] = 24;

   974 	wrap_nc(nc_get_vara_int(input_ncid, mcdate_varid, mc_start, mc_count,

   975 		mcdate));

   976 	wrap_nc(nc_get_vara_int(input_ncid, mcsec_varid, mc_start, mc_count,

   977 		mcsec));

   979 	/* Close the last input file */

   980 	wrap_nc(nc_close(input_ncid));

   982 	/* Store time bounds and compute new time stamps as midpoints */

   983 	for (i = 0; i < 24; i++) {

   984 		tbounds[(i*2)] = bnd1[i];

   985 		tbounds[(i*2+1)] = bnd2[i];

   986 		/* Pick one (integer) day near the midpoint of the bounds and

   987 		 * add in the fraction from the right bounds, so for January

   988 		 * the first hour-of-day value would be 15.0 (plus some year

   989 		 * offset) */

   990 		day1 = floor(bnd1[i]);

   991 		day2 = floor(bnd2[i]);

   992 		times[i] = floor((day1 + day2) / 2.0) + (bnd2[i]-day2);

   993 	}

   995 	return;

   996 }

   998 char **make_cell_methods(size_t nsamples)

   999 {

  1000 	int i;

  1001 	size_t len;

  1002 	char **cell_methods;

  1004 	if (!(cell_methods = (char **)malloc((NUM_TYPES * sizeof(char *))))) {

  1005 		perror("make_cell_methods:cell_methods");

  1006 		exit(2);

  1007 	}

  1008 	for (i = 0; i < NUM_TYPES; i++) {

  1009 		len = strlen(cell_methods_prefix[i]) + 48;

  1010 		if (!(cell_methods[i] = (char *)malloc((len * sizeof(char))))) {

  1011 			perror("make_cell_methods:cell_methods[i]");

  1012 			exit(2);

  1013 		}

  1014 		sprintf(cell_methods[i], "%s (comment: %ld samples)",

  1015 			cell_methods_prefix[i], (long)nsamples);

  1016 	}

  1017 	return cell_methods;

  1018 }

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

  1021 {

  1022 	int i;

  1023 	char *mean_history_gatt, *stddev_history_gatt,

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

  1025 		*stddev_prefix="Statistical standard deviations from history files:",

  1026 		**cell_methods;

  1028 	/*

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

  1030 	 * files.

  1031 	 */

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

  1033 		perror("open_inout:mean_history_gatt");

  1034 		exit(2);

  1035 	}

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

  1037 		perror("open_inout:stddev_history_gatt");

  1038 		exit(2);

  1039 	}

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

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

  1043 	/* Allocate space for time values for each hour of the day */

  1044 	if (!(times = (double *)malloc((24 * sizeof(double))))) {

  1045 		perror("open_inout:times");

  1046 		exit(2);

  1047 	}

  1048 	/* Allocate space for two time bounds values for each hour of the day */

  1049 	if (!(tbounds = (double *)malloc((24 * 2 * sizeof(double))))) {

  1050 		perror("open_inout:tbounds");

  1051 		exit(2);

  1052 	}

  1053 	/* Allocate space for current date and seconds for each hour of day */

  1054 	if (!(mcdate = (int *)malloc((24 * sizeof(int))))) {

  1055 		perror("open_inout:mcdate");

  1056 		exit(2);

  1057 	}

  1058 	if (!(mcsec = (int *)malloc((24 * sizeof(int))))) {

  1059 		perror("open_inout:mcsec");

  1060 		exit(2);

  1061 	}

  1064 	/*

  1065 	 * Explicitly handle "time", "time_bounds" -> "climatology_bounds",

  1066 	 * mcdate, and mcsec metadata variables by constructing their values

  1067 	 * before doing anything else.

  1068 	 */

  1069 	get_time_bounds(first_fname, last_fname, times, tbounds, mcdate, mcsec);

  1070 #ifdef DEBUG

  1071 	for (i = 0; i < 24; i++)

  1072 		fprintf(stderr, "Hour of day=%d, time=%lf, tbounds=%lf,%lf mcdate=%d mcsec=%d\n", i, times[i], tbounds[(i*2)], tbounds[(i*2+1)], mcdate[i], mcsec[i]);

  1073 #endif /* DEBUG */

  1075 	/* Open last input file */

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

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

  1078 	 * and the ID of the unlimited dimension

  1079 	 */

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

  1081 		&input_unlimdimid));

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

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

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

  1085 	get_gatts(input_ncid, input_ngatts, &input_gatt_head);

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

  1087 	get_vars(input_ncid, input_nvars, &input_var_head);

  1088 #ifdef DEBUG

  1089 	varlist_print(input_var_head);

  1090 #endif /* DEBUG */

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

  1092 	/* Create new netCDF files */

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

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

  1095 	/* Define dimensions */

  1096 	mean_ndims = put_dimensions(input_dim_head, input_ndims,

  1097 		input_unlimdimid, nsamples, &idid2mdid, mean_ncid,

  1098 		&mean_dim_head, &mean_unlimdimid);

  1099 	stddev_ndims = put_dimensions(input_dim_head, input_ndims,

  1100 		input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,

  1101 		&stddev_dim_head, &stddev_unlimdimid);

  1102 	/* Define cell_methods[], including number of samples */

  1103 	cell_methods = make_cell_methods(nsamples);

  1104 	/* Define variables and their attributes */

  1105 	mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

  1106 		idid2mdid, mean_ncid, &mean_var_head, cell_methods,

  1107 		MEAN_TYPE);

  1108 	stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,

  1109 		idid2sdid, stddev_ncid, &stddev_var_head, cell_methods,

  1110 		STDDEV_TYPE);

  1111 #ifdef DEBUG

  1112 	printf("Mean variables:\n");

  1113 	varlist_print(mean_var_head);

  1114 	printf("Stddev variables:\n");

  1115 	varlist_print(stddev_var_head);

  1116 #endif /* DEBUG */

  1117 	/* Store global attributes */

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

  1119 	stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,

  1120 		stddev_history_gatt);

  1121 	/* End define mode */

  1122 	wrap_nc(nc_enddef(mean_ncid));

  1123 	wrap_nc(nc_enddef(stddev_ncid));

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

  1125 	copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,

  1126 		mean_var_head, times, tbounds, mcdate, mcsec);

  1127 	copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,

  1128 		stddev_var_head, times, tbounds, mcdate, mcsec);

  1130 	wrap_nc(nc_close(input_ncid));

  1132 	return;

  1133 }

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

  1136 {

  1137 	int i;

  1138 	char name[NC_MAX_NAME+1];

  1139 	size_t len, cnt = 0;

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

  1142 		/* Open input file */

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

  1144 		/*

  1145 		 * Inquire about number of dimensions, variables, global

  1146 		 * attributes and the ID of the unlimited dimension

  1147 		 */

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

  1149 			&input_ngatts, &input_unlimdimid));

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

  1151 		if (strcmp(name, time_name)) {

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

  1153 			exit(4);

  1154 		}

  1155 #ifdef DEBUG

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

  1157 #endif /* DEBUG */

  1158 		wrap_nc(nc_close(input_ncid));

  1159 		cnt += len;

  1160 	}

  1161 	return cnt;

  1162 }

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

  1165 {

  1166 	/*

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

  1168 	 * hour of day.

  1169 	 */

  1170 	int i;

  1171 	size_t **cell_samples;

  1172 	double **means, **stddevs;

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

  1175 		perror("alloc_means_stddevs:*meansp");

  1176 		exit(2);

  1177 	}

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

  1179 		perror("alloc_means_stddevs:*stddevsp");

  1180 		exit(2);

  1181 	}

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

  1183 		perror("alloc_means_stddevs:*cell_samplesp");

  1184 		exit(2);

  1185 	}

  1187 	means = *meansp;

  1188 	stddevs = *stddevsp;

  1189 	cell_samples = *cell_samplesp;

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

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

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

  1194 			exit(2);

  1195 		}

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

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

  1198 			exit(2);

  1199 		}

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

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

  1202 			exit(2);

  1203 		}

  1204 	}

  1206 	return;

  1207 }

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

  1210 {

  1211 	int hr, i, j, pos;

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

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

  1215 #pragma _CRI concurrent

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

  1217 				pos = i*d2+j;

  1218 				means[hr][pos] = FillValue;

  1219 				stddevs[hr][pos] = FillValue;

  1220 				cell_samples[hr][pos] = 0;

  1221 			}

  1222 		}

  1223 	}

  1225 	return;

  1226 }

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

  1229 {

  1230 	int i, j, hr, pos;

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

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

  1234 #pragma _CRI concurrent

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

  1236 				pos = i*d2+j;

  1237 				cell_samples[hr][pos] = 0;

  1238 			}

  1239 		}

  1240 	}

  1242 	return;

  1243 }

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

  1246 {

  1247 	/*

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

  1249 	 * hour of day.

  1250 	 */

  1251 	int i;

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

  1254 		free(means[i]);

  1255 		free(stddevs[i]);

  1256 		free(cell_samples[i]);

  1257 	}

  1259 	free(means);

  1260 	free(stddevs);

  1261 	free(cell_samples);

  1263 	return;

  1264 }

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

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

  1268 {

  1269 	int i, j, hr, pos;

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

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

  1274 #pragma _CRI concurrent

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

  1276 			pos = i*d2+j;

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

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

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

  1280 				else

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

  1282 				++cell_samples[hr][pos];

  1283 			}

  1284 		}

  1285 	}

  1287 	return;

  1288 }

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

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

  1292 {

  1293 	int i, j, hr, pos;

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

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

  1298 #pragma _CRI concurrent

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

  1300 			pos = i*d2+j;

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

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

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

  1304 				else

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

  1306 				++cell_samples[hr][pos];

  1307 			}

  1308 		}

  1309 	}

  1311 	return;

  1312 }

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

  1316 {

  1317 	int i, j, hr, pos;

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

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

  1321 #pragma _CRI concurrent

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

  1323 				pos = i*d2+j;

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

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

  1326 				}

  1327 			}

  1328 		}

  1329 	}

  1331 	return;

  1332 }

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

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

  1336 	size_t **cell_samples)

  1337 {

  1338 	int i, j, hr, pos;

  1339 	double delta;

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

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

  1344 #pragma _CRI concurrent

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

  1346 			pos = i*d2+j;

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

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

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

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

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

  1352 				else

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

  1354 				++cell_samples[hr][pos];

  1355 			}

  1356 		}

  1357 	}

  1359 	return;

  1360 }

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

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

  1364 	size_t **cell_samples)

  1365 {

  1366 	int i, j, hr, pos;

  1367 	double delta;

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

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

  1372 #pragma _CRI concurrent

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

  1374 			pos = i*d2+j;

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

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

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

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

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

  1380 				else

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

  1382 				++cell_samples[hr][pos];

  1383 			}

  1384 		}

  1385 	}

  1387 	return;

  1388 }

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

  1392 {

  1393 	int i, j, hr, pos;

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

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

  1397 #pragma _CRI concurrent

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

  1399 				pos = i*d2+j;

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

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

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

  1403 				}

  1404 			}

  1405 		}

  1406 	}

  1408 	return;

  1409 }

  1411 void read_and_compute(int nifnames, char **input_fnames, size_t d1, size_t d2, char *var_name, nc_type nctype, int level, int ndims, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples)

  1412 {

  1413 	int i, ts;

  1414 	int varid, *mcsec;

  1415 	size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];

  1416 	char name[NC_MAX_NAME+1];

  1417 	void *val;

  1419 	/* Allocate space for timeslice */

  1420 	val = alloc_var(nctype, (d1 * d2));

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

  1423 #ifdef DEBUG

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

  1425 		if (ndims > 3) printf(" to retrieve level %d", level);

  1426 		if (stat_type == MEAN_TYPE)

  1427 			printf(" for computing means\n");

  1428 		else

  1429 			printf(" for computing stddevs\n");

  1430 #endif /* DEBUG */

  1431 		/* Open input file */

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

  1433 		/*

  1434 		 * Inquire about number of dimensions, variables, global

  1435 		 * attributes and the ID of the unlimited dimension

  1436 		 */

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

  1438 			&input_ngatts, &input_unlimdimid));

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

  1440 		if (strcmp(name, time_name)) {

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

  1442 			exit(4);

  1443 		}

  1444 		/* Get seconds of day variable */

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

  1446 		if (!(mcsec = (int *)malloc(len * sizeof(int)))) {

  1447 			perror("read_and_compute:mcsec");

  1448 			exit(2);

  1449 		}

  1450 		wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));

  1451 		/* Get variable ID for requested variable */

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

  1453 		/* Retrieve time slice of desired variable */

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

  1455 			start[0] = ts;

  1456 			count[0] = 1;

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

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

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

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

  1461 			if (ndims == 4) {

  1462 				start[1] = level;

  1463 				count[1] = 1;

  1464 			}

  1465 			switch(nctype) {

  1466 				case NC_FLOAT:

  1467 					wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));

  1468 					if (stat_type == MEAN_TYPE)

  1469 						add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);

  1470 					else

  1471 						add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);

  1472 					break;

  1473 				case NC_DOUBLE:

  1474 					wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));

  1475 					if (stat_type == MEAN_TYPE)

  1476 						add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);

  1477 					else

  1478 						add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);

  1479 					break;

  1480 				default:

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

  1482 					exit(3);

  1483 			}

  1484 		}

  1486 		/* Free mcsec vector */

  1487 		free(mcsec);

  1488 		/* Close input file */

  1489 		wrap_nc(nc_close(input_ncid));

  1490 	}

  1491 	/* Divide sums by number of samples */

  1492 	if (stat_type == MEAN_TYPE)

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

  1494 	else

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

  1497 	/* Free working variable array */

  1498 	free(val);

  1500 	return;

  1501 }

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

  1504 {

  1505 	int i;

  1506 	float *fvar;

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

  1509 		perror("double_to_float:fvar");

  1510 		exit(2);

  1511 	}

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

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

  1516 	return fvar;

  1517 }

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

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

  1521 {

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

  1523 	int i, hr;

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

  1525 	float *fvar = NULL;

  1527 	if (nctype == NC_FLOAT) {

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

  1529 			perror("write_var_hours:fvar");

  1530 			exit(2);

  1531 		}

  1532 	}

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

  1535 		start[0] = hr;

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

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

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

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

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

  1541 		if (ndims == 4) {

  1542 			start[1] = level;

  1543 			count[1] = 1;

  1544 		}

  1545 		switch (nctype) {

  1546 			case NC_FLOAT:

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

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

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

  1550 				break;

  1551 			case NC_DOUBLE:

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

  1553 				break;

  1554 			default:

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

  1556 				exit(3);

  1557 		}

  1558 	}

  1560 	if (nctype == NC_FLOAT)

  1561 		free(fvar);

  1563 	return;

  1564 }

  1566 void compute_stats(int nifnames, char **input_fnames)

  1567 {

  1568 	int j, k, nlevels;

  1569 	size_t d1, d2, **cell_samples;

  1570 	double **means, **stddevs;

  1571 	struct var *in_vnode, *out_vnode;

  1573 	if (input_var_head) {

  1574 		in_vnode = input_var_head;

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

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

  1577 			if (!in_vnode->metadata) {

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

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

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

  1581 					exit(3);

  1582 				}

  1583 				/* Figure out input dimensions */

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

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

  1586 					exit(3);

  1587 				}

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

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

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

  1591 				if (in_vnode->ndims == 3)

  1592 					nlevels = 1;

  1593 				else

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

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

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

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

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

  1599 					in_vnode->name);

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

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

  1602 					/* Read and compute means */

  1603 					read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);

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

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

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

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

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

  1609 					reset_cell_samples(d1, d2, cell_samples);

  1610 					/* Read and compute stddevs using means */

  1611 					read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);

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

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

  1614 					/* Write out stddevs for this variable */

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

  1616 				}

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

  1619 				free_means_stddevs(means, stddevs, cell_samples);

  1620 			}

  1621 			in_vnode = in_vnode->next;

  1622 		}

  1623 	}

  1624 	return;

  1625 }

  1627 void usage(char *program)

  1628 {

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

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

  1631 	exit(4);

  1632 }

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

  1635 {

  1636 	int i, ic, nifnames;

  1637 	size_t len, pos, nsamples;

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

  1640 	ifnames = NULL;

  1641 	mfname = sfname = flist = NULL;

  1642 	nifnames = 0;

  1644 #ifdef DEBUG

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

  1646 	fflush(stdout);

  1647 #endif /* DEBUG */

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

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

  1651 		switch(ic) {

  1652 			case 'm':

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

  1654 					perror("mfname");

  1655 					exit(2);

  1656 				}

  1657 				break;

  1658 			case 's':

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

  1660 					perror("sfname");

  1661 					exit(2);

  1662 				}

  1663 				break;

  1664 		}

  1665 	}

  1667 	if (!mfname) {

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

  1669 		usage(argv[0]);

  1670 	}

  1671 	if (!sfname) {

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

  1673 		usage(argv[0]);

  1674 	}

  1676 	if (optind < argc) {

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

  1678 			perror("ifnames");

  1679 			exit(2);

  1680 		}

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

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

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

  1684 				exit(2);

  1685 			}

  1686 		}

  1687 	}

  1688 	else {

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

  1690 		usage(argv[0]);

  1691 	}

  1694 	/*

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

  1696 	 * attribute on the two outputs files.

  1697 	 */

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

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

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

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

  1702 		perror("flist");

  1703 		exit(2);

  1704 	}

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

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

  1707 #ifdef DEBUG

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

  1709 	fflush(stdout);

  1710 #endif /* DEBUG */

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

  1713 	nsamples = count_samples(nifnames, ifnames);

  1714 #ifdef DEBUG

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

  1716 	fflush(stdout);

  1717 #endif /* DEBUG */

  1719 	/*

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

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

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

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

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

  1725 	 * the statistics are computed.

  1726 	 */

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

  1728 		nsamples);

  1730 	compute_stats(nifnames, ifnames);

  1732 	/* Close the two output files */

  1733 	wrap_nc(nc_close(mean_ncid));

  1734 	wrap_nc(nc_close(stddev_ncid));

  1736 	return 0;

  1737 }

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