1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/h1_summary.c Wed Sep 26 17:16:40 2007 -0400
1.3 @@ -0,0 +1,1525 @@
1.4 +#include <stdio.h>
1.5 +#include <unistd.h>
1.6 +#include <stdlib.h>
1.7 +#include <string.h>
1.8 +#include <math.h>
1.9 +#include "netcdf.h"
1.10 +
1.11 +/*
1.12 + * Loop through one month of hourly history tapes from the Community Land Model
1.13 + * and generate monthly statistics (means and standard deviations) of fields
1.14 + * by hour of day.
1.15 + */
1.16 +
1.17 +#define MEAN_TYPE 0
1.18 +#define STDDEV_TYPE 1
1.19 +
1.20 +#define HOURS_PER_DAY 24
1.21 +#define MIN_PER_HOUR 60
1.22 +#define SEC_PER_MIN 60
1.23 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
1.24 +
1.25 +static char *metadata_vars[] = {
1.26 + "lon",
1.27 + "lat",
1.28 + "lonrof",
1.29 + "latrof",
1.30 + "time",
1.31 + "hour", /* new metadata variable to be added to output files */
1.32 + "levsoi",
1.33 + "levlak",
1.34 + "edgen",
1.35 + "edgee",
1.36 + "edges",
1.37 + "edgew",
1.38 + "longxy",
1.39 + "latixy",
1.40 + "area",
1.41 + "areaupsc",
1.42 + "topo",
1.43 + "topodnsc",
1.44 + "landfrac",
1.45 + "landmask",
1.46 + "pftmask",
1.47 + "indxupsc",
1.48 + "mcdate",
1.49 + "mcsec",
1.50 + "mdcur",
1.51 + "mscur",
1.52 + "nstep",
1.53 + "time_bounds",
1.54 + "date_written",
1.55 + "time_written"
1.56 +};
1.57 +
1.58 +struct text_att {
1.59 + char *name;
1.60 + char *value;
1.61 + struct text_att *next;
1.62 + struct text_att *prev;
1.63 +};
1.64 +
1.65 +struct dim {
1.66 + int dimid;
1.67 + char *name;
1.68 + size_t len;
1.69 + struct dim *next;
1.70 + struct dim *prev;
1.71 +};
1.72 +
1.73 +struct var {
1.74 + int ncvarid;
1.75 + char *name;
1.76 + nc_type nctype;
1.77 + int ndims;
1.78 + int *dimids;
1.79 + int natts;
1.80 + char metadata;
1.81 + char FillFlag;
1.82 + float FillValue;
1.83 + struct var *next;
1.84 + struct var *prev;
1.85 +};
1.86 +
1.87 +static char *time_name = "time";
1.88 +static char *mcsec_name = "mcsec";
1.89 +static char *history_name = "history";
1.90 +static char *nsamples_name = "nsamples";
1.91 +static char *hour_name = "hour", *hour_long_name = "hour of day",
1.92 + *hour_units = "hour";
1.93 +static char *cell_method_name = "cell_method";
1.94 +/* input stuff */
1.95 +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
1.96 +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
1.97 +static struct text_att *input_gatt_head = NULL;
1.98 +static struct dim *input_dim_head = NULL, **input_dim_idx;
1.99 +static struct var *input_var_head = NULL;
1.100 +/* translation stuff */
1.101 +static int *idid2mdid, *idid2sdid; /* dimension IDs */
1.102 +/* output stuff */
1.103 +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
1.104 +static int mean_hour_dimid; /* special notes */
1.105 +static struct dim *mean_dim_head = NULL;
1.106 +static struct var *mean_var_head = NULL;
1.107 +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
1.108 +static int stddev_hour_dimid; /* special notes */
1.109 +static struct dim *stddev_dim_head = NULL;
1.110 +static struct var *stddev_var_head = NULL;
1.111 +
1.112 +char is_metadata(char *name)
1.113 +{
1.114 + int i;
1.115 +
1.116 + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
1.117 +
1.118 + if (i < nmvars)
1.119 + return 1;
1.120 + else
1.121 + return 0;
1.122 +}
1.123 +
1.124 +#if 0
1.125 +struct dim *dimlist_find_by_name(struct dim *head, char *name)
1.126 +{
1.127 + int i;
1.128 + struct dim *p = NULL;
1.129 +
1.130 + if (head) {
1.131 + node = head;
1.132 + for (i = 0 ; i == 0 || node != head; i++) {
1.133 + if (!strcmp(node->name, name))
1.134 + p = node;
1.135 + node = node->next;
1.136 + }
1.137 + return p;
1.138 + }
1.139 + else
1.140 + return NULL;
1.141 +}
1.142 +#endif
1.143 +
1.144 +struct var *varlist_find_by_name(struct var *head, char *name)
1.145 +{
1.146 + int i;
1.147 + struct var *node;
1.148 +
1.149 + if (head) {
1.150 + node = head;
1.151 + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
1.152 + node = node->next;
1.153 + if (i && node == head)
1.154 + return NULL;
1.155 + else
1.156 + return node;
1.157 + }
1.158 + else
1.159 + return NULL;
1.160 +}
1.161 +
1.162 +void gattlist_add_node(struct text_att **headp, char *name, char *value)
1.163 +{
1.164 + struct text_att *head, *node;
1.165 +
1.166 + head = *headp;
1.167 +
1.168 + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
1.169 + perror("gattlist_add_node");
1.170 + exit(2);
1.171 + }
1.172 +
1.173 + if (!(node->name = strdup(name))) {
1.174 + perror("gattlist_add_node");
1.175 + exit(2);
1.176 + }
1.177 + if (!(node->value = strdup(value))) {
1.178 + perror("gattlist_add_node");
1.179 + exit(2);
1.180 + }
1.181 +
1.182 + if (head == NULL) {
1.183 + /* first one */
1.184 + *headp = node;
1.185 + node->next = node;
1.186 + node->prev = node;
1.187 + }
1.188 + else {
1.189 + node->next = head;
1.190 + node->prev = head->prev;
1.191 + /* set this after setting node->prev from here */
1.192 + head->prev = node;
1.193 + /* set this after having set node->prev */
1.194 + node->prev->next = node;
1.195 + }
1.196 +
1.197 + return;
1.198 +}
1.199 +
1.200 +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
1.201 +{
1.202 + struct dim *head, *node;
1.203 +
1.204 + head = *headp;
1.205 +
1.206 + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
1.207 + perror("dimlist_add_node");
1.208 + exit(2);
1.209 + }
1.210 +
1.211 + node->dimid = dimid;
1.212 + if (!(node->name = strdup(name))) {
1.213 + perror("dimlist_add_node");
1.214 + exit(2);
1.215 + }
1.216 + node->len = len;
1.217 +
1.218 + if (head == NULL) {
1.219 + /* first one */
1.220 + *headp = node;
1.221 + node->next = node;
1.222 + node->prev = node;
1.223 + }
1.224 + else {
1.225 + node->next = head;
1.226 + node->prev = head->prev;
1.227 + /* set this after setting node->prev from here */
1.228 + head->prev = node;
1.229 + /* set this after having set node->prev */
1.230 + node->prev->next = node;
1.231 + }
1.232 +
1.233 + return node;
1.234 +}
1.235 +
1.236 +void varlist_add_node(struct var **headp, int ncvarid, char *name,
1.237 + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
1.238 + float FillValue)
1.239 +{
1.240 + int i;
1.241 + struct var *head, *node;
1.242 +
1.243 + head = *headp;
1.244 +
1.245 + if (!(node = (struct var *)malloc(sizeof(struct var)))) {
1.246 + perror("varlist_add_node");
1.247 + exit(3);
1.248 + }
1.249 +
1.250 + node->ncvarid = ncvarid;
1.251 + if (!(node->name = strdup(name))) {
1.252 + perror("varlist_add_node");
1.253 + exit(3);
1.254 + }
1.255 + node->nctype = nctype;
1.256 + node->ndims = ndims;
1.257 + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
1.258 + perror("varlist_add_node");
1.259 + exit(3);
1.260 + }
1.261 + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
1.262 + node->natts = natts;
1.263 + node->metadata = is_metadata(name);
1.264 + node->FillFlag = FillFlag;
1.265 + node->FillValue = FillValue;
1.266 +
1.267 + if (head == NULL) {
1.268 + /* first one */
1.269 + *headp = node;
1.270 + node->next = node;
1.271 + node->prev = node;
1.272 + }
1.273 + else {
1.274 + node->next = head;
1.275 + node->prev = head->prev;
1.276 + /* set this after setting node->prev from here */
1.277 + head->prev = node;
1.278 + /* set this after having set node->prev */
1.279 + node->prev->next = node;
1.280 + }
1.281 +
1.282 + return;
1.283 +}
1.284 +
1.285 +void varlist_print(struct var *headp)
1.286 +{
1.287 + int i, j;
1.288 + struct var *node;
1.289 +
1.290 + printf("Beginning of Variable List\n");
1.291 +
1.292 + if (headp) {
1.293 + node = headp;
1.294 + for (j = 0; j == 0 || node != headp; j++) {
1.295 + printf("Variable %d (ptr=%ld):\n", j, (long)node);
1.296 + printf("\tncvarid = %d\n", node->ncvarid);
1.297 + printf("\tname = %s\n", node->name);
1.298 + printf("\tnctype = %d\n", node->nctype);
1.299 + printf("\tndims = %d\n", node->ndims);
1.300 + printf("\tdimids =");
1.301 + for (i = 0; i < node->ndims; i++)
1.302 + printf(" %d", node->dimids[i]);
1.303 + printf("\n");
1.304 + printf("\tnatts = %d\n", node->natts);
1.305 + printf("\tmetadata = %d\n", (int)node->metadata);
1.306 + printf("\tnext = %ld\n", (long)node->next);
1.307 + printf("\tprev = %ld\n", (long)node->prev);
1.308 + node = node->next;
1.309 + }
1.310 + }
1.311 + else {
1.312 + printf("\tList undefined\n");
1.313 + }
1.314 +
1.315 + printf("End of Variable List\n");
1.316 +
1.317 + return;
1.318 +}
1.319 +
1.320 +void wrap_nc(int status)
1.321 +{
1.322 + if (status != NC_NOERR) {
1.323 + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
1.324 + exit(1);
1.325 + }
1.326 +
1.327 + return;
1.328 +}
1.329 +
1.330 +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
1.331 +{
1.332 + int i;
1.333 + char name[NC_MAX_NAME+1];
1.334 + size_t len;
1.335 + struct dim **in_dim_idx;
1.336 +
1.337 + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
1.338 + perror("*in_dim_idxp");
1.339 + exit(3);
1.340 + }
1.341 + in_dim_idx = *in_dim_idxp;
1.342 +
1.343 + for (i = 0; i < ndims; i++) {
1.344 + wrap_nc(nc_inq_dim(ncid, i, name, &len));
1.345 + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
1.346 + }
1.347 +
1.348 + return;
1.349 +}
1.350 +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
1.351 +{
1.352 + int i;
1.353 + char name[NC_MAX_NAME+1], *value;
1.354 + nc_type xtype;
1.355 + size_t len;
1.356 +
1.357 + for (i = 0; i < ngatts; i++) {
1.358 + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
1.359 + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
1.360 + if (xtype != NC_CHAR) {
1.361 + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
1.362 + exit(2);
1.363 + }
1.364 + if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
1.365 + perror("get_gatts: value");
1.366 + exit(3);
1.367 + }
1.368 + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
1.369 + value[(len+1-1)] = '\0';
1.370 + gattlist_add_node(gatt_headp, name, value);
1.371 + free(value); /* because gattlist_add_node() duplicates it */
1.372 + }
1.373 +
1.374 + return;
1.375 +}
1.376 +
1.377 +void get_vars(int ncid, int nvars, struct var **var_headp)
1.378 +{
1.379 + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
1.380 + size_t len;
1.381 + float FillValue;
1.382 + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
1.383 + nc_type xtype, atype;
1.384 +
1.385 + for (i = 0; i < nvars; i++) {
1.386 + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
1.387 + &natts));
1.388 + /* Check for _FillValue */
1.389 + FillFlag = 0;
1.390 + FillValue = 0.;
1.391 + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
1.392 + == NC_NOERR) {
1.393 + if (atype == NC_FLOAT && len) {
1.394 + wrap_nc(nc_get_att_float(ncid, i,
1.395 + fill_att_name, &FillValue));
1.396 + FillFlag = 1;
1.397 + }
1.398 + }
1.399 + varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
1.400 + natts, FillFlag, FillValue);
1.401 + }
1.402 +
1.403 + return;
1.404 +}
1.405 +
1.406 +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
1.407 + size_t nsamples, int **in2outp, int out_ncid,
1.408 + struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)
1.409 +{
1.410 + /*
1.411 + * Define dimensions on new files and build translation tables between
1.412 + * dimension IDs.
1.413 + */
1.414 + int j, dimid, ndims, *in2out;
1.415 + size_t len;
1.416 + struct dim *dnode;
1.417 +
1.418 + ndims = 0;
1.419 +
1.420 + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
1.421 + perror("put_dimensions: in2outp");
1.422 + exit(3);
1.423 + }
1.424 + in2out = *in2outp;
1.425 +
1.426 + if (in_dim_head) {
1.427 + dnode = in_dim_head;
1.428 + for (j = 0; j == 0 || dnode != in_dim_head; j++) {
1.429 + if (dnode->dimid == in_unlimdimid)
1.430 + len = NC_UNLIMITED;
1.431 + else
1.432 + len = dnode->len;
1.433 + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
1.434 + dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
1.435 + ++ndims;
1.436 + in2out[dnode->dimid] = dimid;
1.437 + if (dnode->dimid == in_unlimdimid)
1.438 + *out_unlimdimidp = dimid;
1.439 + /*
1.440 + * Just after the "time" dimension, add the new
1.441 + * "nsamples" and "hour" dimensions.
1.442 + */
1.443 + if (!strcmp(dnode->name, time_name)) {
1.444 + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
1.445 + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
1.446 + ++ndims;
1.447 +
1.448 + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
1.449 + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
1.450 + ++ndims;
1.451 + /* Track hour dimid for out files */
1.452 + *out_hour_dimidp = dimid;
1.453 + }
1.454 +
1.455 + dnode = dnode->next;
1.456 + }
1.457 + }
1.458 + else {
1.459 + fprintf(stderr, "WARNING: No dimensions defined!\n");
1.460 + }
1.461 +
1.462 + return ndims;
1.463 +}
1.464 +
1.465 +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
1.466 +{
1.467 + /*
1.468 + * Write out global attributes matching those of the input file.
1.469 + * Change history attribute to the string passed in as an argument.
1.470 + */
1.471 + int j, hflag = 0, ngatts;
1.472 + char *value;
1.473 + struct text_att *anode;
1.474 +
1.475 + ngatts = 0;
1.476 +
1.477 + if (in_gatt_head) {
1.478 + anode = in_gatt_head;
1.479 + for (j = 0; j == 0 || anode != in_gatt_head; j++) {
1.480 + if (!strcmp(anode->name, history_name)) {
1.481 + value = out_history;
1.482 + ++hflag;
1.483 + }
1.484 + else
1.485 + value = anode->value;
1.486 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
1.487 + ++ngatts;
1.488 + anode = anode->next;
1.489 + }
1.490 + /* If no history attribute on input, add one to the output */
1.491 + if (!hflag) {
1.492 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
1.493 + ++ngatts;
1.494 + }
1.495 + }
1.496 + else {
1.497 + fprintf(stderr, "WARNING: No global attributes defined!\n");
1.498 + }
1.499 +
1.500 + return ngatts;
1.501 +}
1.502 +
1.503 +int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)
1.504 +{
1.505 + /*
1.506 + * Translate between input dimension IDs and those for the output file.
1.507 + * For normal time-based variables, add a new dimension for hour of
1.508 + * day. For metadata variables, do not add new dimensions, even if
1.509 + * it is time-based. Return (possibly new) number of dimensions.
1.510 + */
1.511 + int i, ndims;
1.512 +
1.513 + for (i = ndims = 0; i < in_ndims; i++, ndims++) {
1.514 + out_dimids[ndims] = in2out[in_dimids[i]];
1.515 + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
1.516 + && !metadata) {
1.517 + ndims++;
1.518 + out_dimids[ndims] = hour_dimid;
1.519 + }
1.520 + }
1.521 +
1.522 + return ndims;
1.523 +}
1.524 +
1.525 +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
1.526 + int in_varid, int out_ncid, int out_varid)
1.527 +{
1.528 + /*
1.529 + * Copy the attributes of the input variable to those of the output
1.530 + * variable. If the variable is not a metadata variable, ensure that
1.531 + * the cell_method attribute is set either to "time: mean" or
1.532 + * "time: standard_deviation", depending on the output file type.
1.533 + */
1.534 +
1.535 + int i, natts;
1.536 + char name[NC_MAX_NAME+1], cmflag = 0;
1.537 + char *cell_methods[] = { "time: mean", "time: standard_deviation" };
1.538 +
1.539 + for (i = natts = 0; i < in_natts; i++, natts++) {
1.540 + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
1.541 + if (!strcmp(name, cell_method_name) && !metadata) {
1.542 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
1.543 + cmflag = 1;
1.544 + }
1.545 + else
1.546 + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
1.547 + }
1.548 + /*
1.549 + * If no cell_method attribute was specified for a non-metadata
1.550 + * variable on the input file, add the appropriate cell_method anyway
1.551 + * on the output file.
1.552 + */
1.553 + if (!cmflag && !metadata) {
1.554 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
1.555 + natts++;
1.556 + }
1.557 +
1.558 + return natts;
1.559 +}
1.560 +
1.561 +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
1.562 + int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp,
1.563 + int stat_type)
1.564 +{
1.565 + /*
1.566 + * Define variables on output file and copy attributes from input file.
1.567 + * Return number of variables defined.
1.568 + */
1.569 + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
1.570 + struct var *vnode;
1.571 +
1.572 + nvars = 0;
1.573 +
1.574 + if (in_var_head) {
1.575 + vnode = in_var_head;
1.576 + /*
1.577 + * March through input variables creating (mostly) the same
1.578 + * variables on the output file.
1.579 + */
1.580 + for (j = 0; j == 0 || vnode != in_var_head; j++) {
1.581 + ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
1.582 + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
1.583 + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
1.584 + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
1.585 + ++nvars;
1.586 + /*
1.587 + * Just after the "time" variable, add the new "hour"
1.588 + * variable for hour of day.
1.589 + */
1.590 + if (!strcmp(vnode->name, time_name)) {
1.591 + ndims = 1;
1.592 + dimids[0] = hour_dimid;
1.593 + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
1.594 + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
1.595 + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
1.596 + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
1.597 + ++nvars;
1.598 + }
1.599 +
1.600 + vnode = vnode->next;
1.601 + }
1.602 + }
1.603 + else {
1.604 + fprintf(stderr, "WARNING: No variables defined!\n");
1.605 + }
1.606 +
1.607 + return nvars;
1.608 +}
1.609 +
1.610 +void *alloc_var(nc_type nctype, size_t len)
1.611 +{
1.612 + void *val;
1.613 +
1.614 + switch (nctype) {
1.615 + case NC_FLOAT:
1.616 + if (!(val = (float *)malloc((len) * sizeof(float)))) {
1.617 + perror("alloc_var: val");
1.618 + exit(3);
1.619 + }
1.620 + break;
1.621 + case NC_INT:
1.622 + if (!(val = (int *)malloc((len) * sizeof(int)))) {
1.623 + perror("alloc_var: val");
1.624 + exit(3);
1.625 + }
1.626 + break;
1.627 + case NC_DOUBLE:
1.628 + if (!(val = (double *)malloc((len) * sizeof(double)))) {
1.629 + perror("alloc_var: val");
1.630 + exit(3);
1.631 + }
1.632 + break;
1.633 + case NC_CHAR:
1.634 + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
1.635 + perror("alloc_var: val");
1.636 + exit(3);
1.637 + }
1.638 + break;
1.639 + default:
1.640 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.641 + exit(3);
1.642 + }
1.643 +
1.644 + return val;
1.645 +}
1.646 +
1.647 +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
1.648 +{
1.649 + int i;
1.650 + size_t len = (size_t)1;
1.651 + void *val;
1.652 +
1.653 + /* Figure out the total size */
1.654 + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
1.655 +
1.656 + val = alloc_var(nctype, len);
1.657 +
1.658 + switch (nctype) {
1.659 + case NC_FLOAT:
1.660 + wrap_nc(nc_get_var_float(ncid, varid, val));
1.661 + break;
1.662 + case NC_INT:
1.663 + wrap_nc(nc_get_var_int(ncid, varid, val));
1.664 + break;
1.665 + case NC_DOUBLE:
1.666 + wrap_nc(nc_get_var_double(ncid, varid, val));
1.667 + break;
1.668 + case NC_CHAR:
1.669 + wrap_nc(nc_get_var_text(ncid, varid, val));
1.670 + break;
1.671 + default:
1.672 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.673 + exit(3);
1.674 + }
1.675 +
1.676 + return val;
1.677 +}
1.678 +
1.679 +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
1.680 +{
1.681 + int i;
1.682 + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1.683 + void *val;
1.684 +
1.685 + /* Make sure time is really the first dimension */
1.686 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
1.687 + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
1.688 + exit(3);
1.689 + }
1.690 + /*
1.691 + * Figure out the total size for the slice (start at second dimension)
1.692 + * and build start/count vectors.
1.693 + */
1.694 + start[0] = tslice;
1.695 + count[0] = 1;
1.696 + for (i = 1; i < ndims; i++) {
1.697 + start[i] = 0;
1.698 + count[i] = (dim_idx[dimids[i]])->len;
1.699 + len *= (dim_idx[dimids[i]])->len;
1.700 + }
1.701 +
1.702 + val = alloc_var(nctype, len);
1.703 +
1.704 + switch (nctype) {
1.705 + case NC_FLOAT:
1.706 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
1.707 + break;
1.708 + case NC_INT:
1.709 + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
1.710 + break;
1.711 + case NC_DOUBLE:
1.712 + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
1.713 + break;
1.714 + case NC_CHAR:
1.715 + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
1.716 + break;
1.717 + default:
1.718 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.719 + exit(3);
1.720 + }
1.721 +
1.722 + return val;
1.723 +}
1.724 +
1.725 +void write_var(int ncid, int varid, nc_type nctype, void *val)
1.726 +{
1.727 + switch (nctype) {
1.728 + case NC_FLOAT:
1.729 + wrap_nc(nc_put_var_float(ncid, varid, val));
1.730 + break;
1.731 + case NC_INT:
1.732 + wrap_nc(nc_put_var_int(ncid, varid, val));
1.733 + break;
1.734 + case NC_DOUBLE:
1.735 + wrap_nc(nc_put_var_double(ncid, varid, val));
1.736 + break;
1.737 + case NC_CHAR:
1.738 + wrap_nc(nc_put_var_text(ncid, varid, val));
1.739 + break;
1.740 + default:
1.741 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.742 + exit(3);
1.743 + }
1.744 +
1.745 + return;
1.746 +}
1.747 +
1.748 +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
1.749 +{
1.750 + int i;
1.751 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1.752 +
1.753 + /* Make sure time is really the first dimension */
1.754 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
1.755 + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
1.756 + exit(3);
1.757 + }
1.758 +
1.759 + /* Build start/count vectors */
1.760 + start[0] = tslice;
1.761 + count[0] = 1;
1.762 + for (i = 1; i < ndims; i++) {
1.763 + start[i] = 0;
1.764 + count[i] = (dim_idx[dimids[i]])->len;
1.765 + }
1.766 +
1.767 + switch (nctype) {
1.768 + case NC_FLOAT:
1.769 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
1.770 + break;
1.771 + case NC_INT:
1.772 + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
1.773 + break;
1.774 + case NC_DOUBLE:
1.775 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
1.776 + break;
1.777 + case NC_CHAR:
1.778 + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
1.779 + break;
1.780 + default:
1.781 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.782 + exit(3);
1.783 + }
1.784 +
1.785 + return;
1.786 +}
1.787 +
1.788 +void copy_metadata(int in_ncid, struct var *in_var_head,
1.789 + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head)
1.790 +{
1.791 + int i, j;
1.792 + float hr[HOURS_PER_DAY];
1.793 + struct var *in_vnode, *out_vnode;
1.794 + void *val;
1.795 +
1.796 + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
1.797 +
1.798 + if (in_var_head) {
1.799 + in_vnode = in_var_head;
1.800 + /*
1.801 + * March through input variables to stuff metadata values into
1.802 + * the new output files. NOTE: Time-based metadata variables
1.803 + * should contain only the last (time-slice) value (from all
1.804 + * input files).
1.805 + */
1.806 + for (j = 0; j == 0 || in_vnode != in_var_head; j++) {
1.807 + if (in_vnode->metadata) {
1.808 + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
1.809 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
1.810 + /* time is not the first dimension */
1.811 +#ifdef DEBUG
1.812 + printf("Copying metadata variable %s\n",
1.813 + in_vnode->name);
1.814 +#endif /* DEBUG */
1.815 + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
1.816 + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
1.817 + }
1.818 + else {
1.819 + /* time is the first dimension */
1.820 +#ifdef DEBUG
1.821 + printf("Copying last value of \
1.822 +time-based metadata variable %s\n", in_vnode->name);
1.823 +#endif /* DEBUG */
1.824 + 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));
1.825 + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
1.826 + }
1.827 + free(val);
1.828 + /*
1.829 + * Just after the "time" variable, write out
1.830 + * the "hour" variable values.
1.831 + */
1.832 + if (!strcmp(in_vnode->name, time_name)) {
1.833 + out_vnode = varlist_find_by_name(out_var_head, hour_name);
1.834 + write_var(out_ncid, out_vnode->ncvarid,
1.835 + out_vnode->nctype, hr);
1.836 + }
1.837 + }
1.838 + else {
1.839 +#ifdef DEBUG
1.840 + printf("Skipping variable %s\n",
1.841 + in_vnode->name);
1.842 +#endif /* DEBUG */
1.843 + }
1.844 + in_vnode = in_vnode->next;
1.845 + }
1.846 + }
1.847 +
1.848 + return;
1.849 +}
1.850 +
1.851 +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
1.852 +{
1.853 + char *mean_history_gatt, *stddev_history_gatt,
1.854 + *mean_prefix="Statistical means from history files:",
1.855 + *stddev_prefix="Statistical standard deviations from history files:";
1.856 +
1.857 + /*
1.858 + * Construct strings for history global attributes for the two output
1.859 + * files.
1.860 + */
1.861 + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
1.862 + perror("open_inout:mean_history_gatt");
1.863 + exit(2);
1.864 + }
1.865 + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
1.866 + perror("open_inout:stddev_history_gatt");
1.867 + exit(2);
1.868 + }
1.869 + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
1.870 + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
1.871 +
1.872 + /* Open input file */
1.873 + wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid));
1.874 + /* Inquire about number of dimensions, variables, global attributes
1.875 + * and the ID of the unlimited dimension
1.876 + */
1.877 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
1.878 + &input_unlimdimid));
1.879 + /* Grab dimension IDs and lengths from input file */
1.880 + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
1.881 + /* Grab desired global attributes from input file */
1.882 + get_gatts(input_ncid, input_ngatts, &input_gatt_head);
1.883 + /* Grab variable IDs and attributes from input file */
1.884 + get_vars(input_ncid, input_nvars, &input_var_head);
1.885 +#ifdef DEBUG
1.886 + varlist_print(input_var_head);
1.887 +#endif /* DEBUG */
1.888 + /* Create netCDF files for monthly means and standard deviations */
1.889 + /* Create new netCDF files */
1.890 + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
1.891 + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
1.892 + /* Define dimensions */
1.893 + mean_ndims = put_dimensions(input_dim_head, input_ndims,
1.894 + input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
1.895 + &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
1.896 + stddev_ndims = put_dimensions(input_dim_head, input_ndims,
1.897 + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
1.898 + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
1.899 + /* Define variables and their attributes */
1.900 + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
1.901 + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
1.902 + MEAN_TYPE);
1.903 + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
1.904 + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
1.905 + STDDEV_TYPE);
1.906 + /* Store global attributes */
1.907 + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
1.908 + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
1.909 + stddev_history_gatt);
1.910 + /* End define mode */
1.911 + wrap_nc(nc_enddef(mean_ncid));
1.912 + wrap_nc(nc_enddef(stddev_ncid));
1.913 + /* Write out metdata variables that do not depend on "time" */
1.914 + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
1.915 + mean_var_head);
1.916 + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
1.917 + stddev_var_head);
1.918 +
1.919 + wrap_nc(nc_close(input_ncid));
1.920 +
1.921 + return;
1.922 +}
1.923 +
1.924 +size_t count_samples(int nifnames, char **input_fnames)
1.925 +{
1.926 + int i;
1.927 + char name[NC_MAX_NAME+1];
1.928 + size_t len, cnt = 0;
1.929 +
1.930 + for (i = 0; i < nifnames; i++) {
1.931 + /* Open input file */
1.932 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1.933 + /*
1.934 + * Inquire about number of dimensions, variables, global
1.935 + * attributes and the ID of the unlimited dimension
1.936 + */
1.937 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1.938 + &input_ngatts, &input_unlimdimid));
1.939 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1.940 + if (strcmp(name, time_name)) {
1.941 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1.942 + exit(4);
1.943 + }
1.944 +#ifdef DEBUG
1.945 + printf("%ld samples in %s\n", (long)len, input_fnames[i]);
1.946 +#endif /* DEBUG */
1.947 + wrap_nc(nc_close(input_ncid));
1.948 + cnt += len;
1.949 + }
1.950 + return cnt;
1.951 +}
1.952 +
1.953 +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
1.954 +{
1.955 + /*
1.956 + * Allocate space for arrays of means and standard deviations by
1.957 + * hour of day.
1.958 + */
1.959 + int i;
1.960 + size_t **cell_samples;
1.961 + double **means, **stddevs;
1.962 +
1.963 + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1.964 + perror("alloc_means_stddevs:*meansp");
1.965 + exit(2);
1.966 + }
1.967 + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
1.968 + perror("alloc_means_stddevs:*stddevsp");
1.969 + exit(2);
1.970 + }
1.971 + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
1.972 + perror("alloc_means_stddevs:*cell_samplesp");
1.973 + exit(2);
1.974 + }
1.975 +
1.976 + means = *meansp;
1.977 + stddevs = *stddevsp;
1.978 + cell_samples = *cell_samplesp;
1.979 +
1.980 + for (i = 0; i < HOURS_PER_DAY; i++) {
1.981 + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1.982 + perror("alloc_means_stddevs:means[i]");
1.983 + exit(2);
1.984 + }
1.985 + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
1.986 + perror("alloc_means_stddevs:stddevs[i]");
1.987 + exit(2);
1.988 + }
1.989 + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
1.990 + perror("alloc_means_stddevs:cell_samples[i]");
1.991 + exit(2);
1.992 + }
1.993 + }
1.994 +
1.995 + return;
1.996 +}
1.997 +
1.998 +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
1.999 +{
1.1000 + int hr, i, j, pos;
1.1001 +
1.1002 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1.1003 + for (i = 0; i < d1; i++) {
1.1004 +#pragma _CRI concurrent
1.1005 + for (j = 0; j < d2; j++) {
1.1006 + pos = i*d2+j;
1.1007 + means[hr][pos] = FillValue;
1.1008 + stddevs[hr][pos] = FillValue;
1.1009 + cell_samples[hr][pos] = 0;
1.1010 + }
1.1011 + }
1.1012 + }
1.1013 +
1.1014 + return;
1.1015 +}
1.1016 +
1.1017 +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
1.1018 +{
1.1019 + int i, j, hr, pos;
1.1020 +
1.1021 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1.1022 + for (i = 0; i < d1; i++) {
1.1023 +#pragma _CRI concurrent
1.1024 + for (j = 0; j < d2; j++) {
1.1025 + pos = i*d2+j;
1.1026 + cell_samples[hr][pos] = 0;
1.1027 + }
1.1028 + }
1.1029 + }
1.1030 +
1.1031 + return;
1.1032 +}
1.1033 +
1.1034 +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
1.1035 +{
1.1036 + /*
1.1037 + * Free space from arrays of means and standard deviations by
1.1038 + * hour of day.
1.1039 + */
1.1040 + int i;
1.1041 +
1.1042 + for (i = 0; i < HOURS_PER_DAY; i++) {
1.1043 + free(means[i]);
1.1044 + free(stddevs[i]);
1.1045 + free(cell_samples[i]);
1.1046 + }
1.1047 +
1.1048 + free(means);
1.1049 + free(stddevs);
1.1050 + free(cell_samples);
1.1051 +
1.1052 + return;
1.1053 +}
1.1054 +
1.1055 +void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
1.1056 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
1.1057 +{
1.1058 + int i, j, hr, pos;
1.1059 +
1.1060 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1.1061 +
1.1062 + for (i = 0; i < d1; i++) {
1.1063 +#pragma _CRI concurrent
1.1064 + for (j = 0; j < d2; j++) {
1.1065 + pos = i*d2+j;
1.1066 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1.1067 + if (cell_samples[hr][pos] == 0)
1.1068 + means[hr][pos] = (double)val[pos];
1.1069 + else
1.1070 + means[hr][pos] += (double)val[pos];
1.1071 + ++cell_samples[hr][pos];
1.1072 + }
1.1073 + }
1.1074 + }
1.1075 +
1.1076 + return;
1.1077 +}
1.1078 +
1.1079 +void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
1.1080 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
1.1081 +{
1.1082 + int i, j, hr, pos;
1.1083 +
1.1084 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1.1085 +
1.1086 + for (i = 0; i < d1; i++) {
1.1087 +#pragma _CRI concurrent
1.1088 + for (j = 0; j < d2; j++) {
1.1089 + pos = i*d2+j;
1.1090 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
1.1091 + if (cell_samples[hr][pos] == 0)
1.1092 + means[hr][pos] = val[pos];
1.1093 + else
1.1094 + means[hr][pos] += val[pos];
1.1095 + ++cell_samples[hr][pos];
1.1096 + }
1.1097 + }
1.1098 + }
1.1099 +
1.1100 + return;
1.1101 +}
1.1102 +
1.1103 +
1.1104 +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
1.1105 +{
1.1106 + int i, j, hr, pos;
1.1107 +
1.1108 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1.1109 + for (i = 0; i < d1; i++) {
1.1110 +#pragma _CRI concurrent
1.1111 + for (j = 0; j < d2; j++) {
1.1112 + pos = i*d2+j;
1.1113 + if (cell_samples[hr][pos] != 0) {
1.1114 + means[hr][pos] /= (double)cell_samples[hr][pos];
1.1115 + }
1.1116 + }
1.1117 + }
1.1118 + }
1.1119 +
1.1120 + return;
1.1121 +}
1.1122 +
1.1123 +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
1.1124 + char FillFlag, float FillValue, double **means, double **stddevs,
1.1125 + size_t **cell_samples)
1.1126 +{
1.1127 + int i, j, hr, pos;
1.1128 + double delta;
1.1129 +
1.1130 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1.1131 +
1.1132 + for (i = 0; i < d1; i++) {
1.1133 +#pragma _CRI concurrent
1.1134 + for (j = 0; j < d2; j++) {
1.1135 + pos = i*d2+j;
1.1136 + if (!FillFlag || (FillFlag && val[pos] != FillValue
1.1137 + && means[hr][pos] != FillValue)) {
1.1138 + delta = means[hr][pos] - (double)val[pos];
1.1139 + if (cell_samples[hr][pos] == 0)
1.1140 + stddevs[hr][pos] = delta * delta;
1.1141 + else
1.1142 + stddevs[hr][pos] += delta * delta;
1.1143 + ++cell_samples[hr][pos];
1.1144 + }
1.1145 + }
1.1146 + }
1.1147 +
1.1148 + return;
1.1149 +}
1.1150 +
1.1151 +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
1.1152 + char FillFlag, float FillValue, double **means, double **stddevs,
1.1153 + size_t **cell_samples)
1.1154 +{
1.1155 + int i, j, hr, pos;
1.1156 + double delta;
1.1157 +
1.1158 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
1.1159 +
1.1160 + for (i = 0; i < d1; i++) {
1.1161 +#pragma _CRI concurrent
1.1162 + for (j = 0; j < d2; j++) {
1.1163 + pos = i*d2+j;
1.1164 + if (!FillFlag || (FillFlag && val[pos] != FillValue
1.1165 + && means[hr][pos] != FillValue)) {
1.1166 + delta = means[hr][pos] - val[pos];
1.1167 + if (cell_samples[hr][pos] == 0)
1.1168 + stddevs[hr][pos] = delta * delta;
1.1169 + else
1.1170 + stddevs[hr][pos] += delta * delta;
1.1171 + ++cell_samples[hr][pos];
1.1172 + }
1.1173 + }
1.1174 + }
1.1175 +
1.1176 + return;
1.1177 +}
1.1178 +
1.1179 +
1.1180 +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
1.1181 +{
1.1182 + int i, j, hr, pos;
1.1183 +
1.1184 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1.1185 + for (i = 0; i < d1; i++) {
1.1186 +#pragma _CRI concurrent
1.1187 + for (j = 0; j < d2; j++) {
1.1188 + pos = i*d2+j;
1.1189 + if (cell_samples[hr][pos] != 0) {
1.1190 + stddevs[hr][pos] /= (double)cell_samples[hr][pos];
1.1191 + stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
1.1192 + }
1.1193 + }
1.1194 + }
1.1195 + }
1.1196 +
1.1197 + return;
1.1198 +}
1.1199 +
1.1200 +
1.1201 +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)
1.1202 +{
1.1203 + int i, ts;
1.1204 + int varid, *mcsec;
1.1205 + size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1.1206 + char name[NC_MAX_NAME+1];
1.1207 + void *val;
1.1208 +
1.1209 + /* Allocate space for timeslice */
1.1210 + val = alloc_var(nctype, (d1 * d2));
1.1211 +
1.1212 + for (i = 0; i < nifnames; i++) {
1.1213 +#ifdef DEBUG
1.1214 + printf("\tOpening %s", input_fnames[i]);
1.1215 + if (ndims > 3) printf(" to retrieve level %d", level);
1.1216 + if (stat_type == MEAN_TYPE)
1.1217 + printf(" for computing means\n");
1.1218 + else
1.1219 + printf(" for computing stddevs\n");
1.1220 +#endif /* DEBUG */
1.1221 + /* Open input file */
1.1222 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
1.1223 + /*
1.1224 + * Inquire about number of dimensions, variables, global
1.1225 + * attributes and the ID of the unlimited dimension
1.1226 + */
1.1227 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
1.1228 + &input_ngatts, &input_unlimdimid));
1.1229 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
1.1230 + if (strcmp(name, time_name)) {
1.1231 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
1.1232 + exit(4);
1.1233 + }
1.1234 + /* Get seconds of day variable */
1.1235 + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
1.1236 + if (!(mcsec = (int *)malloc(len * sizeof(int)))) {
1.1237 + perror("read_and_compute:mcsec");
1.1238 + exit(2);
1.1239 + }
1.1240 + wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));
1.1241 + /* Get variable ID for requested variable */
1.1242 + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
1.1243 + /* Retrieve time slice of desired variable */
1.1244 + for (ts = 0; ts < len; ts++) {
1.1245 + start[0] = ts;
1.1246 + count[0] = 1;
1.1247 + start[(ndims-2)] = 0;
1.1248 + count[(ndims-2)] = d1;
1.1249 + start[(ndims-1)] = 0;
1.1250 + count[(ndims-1)] = d2;
1.1251 + if (ndims == 4) {
1.1252 + start[1] = level;
1.1253 + count[1] = 1;
1.1254 + }
1.1255 + switch(nctype) {
1.1256 + case NC_FLOAT:
1.1257 + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));
1.1258 + if (stat_type == MEAN_TYPE)
1.1259 + add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1.1260 + else
1.1261 + add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1.1262 + break;
1.1263 + case NC_DOUBLE:
1.1264 + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));
1.1265 + if (stat_type == MEAN_TYPE)
1.1266 + add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
1.1267 + else
1.1268 + add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
1.1269 + break;
1.1270 + default:
1.1271 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.1272 + exit(3);
1.1273 + }
1.1274 + }
1.1275 +
1.1276 + /* Free mcsec vector */
1.1277 + free(mcsec);
1.1278 + /* Close input file */
1.1279 + wrap_nc(nc_close(input_ncid));
1.1280 + }
1.1281 + /* Divide sums by number of samples */
1.1282 + if (stat_type == MEAN_TYPE)
1.1283 + divide_means(d1, d2, means, cell_samples);
1.1284 + else
1.1285 + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
1.1286 +
1.1287 + /* Free working variable array */
1.1288 + free(val);
1.1289 +
1.1290 + return;
1.1291 +}
1.1292 +
1.1293 +float *double_to_float(size_t len, double *dvar)
1.1294 +{
1.1295 + int i;
1.1296 + float *fvar;
1.1297 +
1.1298 + if (!(fvar = (float *)malloc(len * sizeof(float)))) {
1.1299 + perror("double_to_float:fvar");
1.1300 + exit(2);
1.1301 + }
1.1302 +
1.1303 + for (i = 0; i < len; i++)
1.1304 + fvar[i] = (float)dvar[i];
1.1305 +
1.1306 + return fvar;
1.1307 +}
1.1308 +
1.1309 +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
1.1310 + int level, int ndims, double **var_hours)
1.1311 +{
1.1312 + /* Output dimensions should be one larger than input dimensions */
1.1313 + int i, hr;
1.1314 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
1.1315 + float *fvar = NULL;
1.1316 +
1.1317 + if (nctype == NC_FLOAT) {
1.1318 + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
1.1319 + perror("write_var_hours:fvar");
1.1320 + exit(2);
1.1321 + }
1.1322 + }
1.1323 +
1.1324 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
1.1325 + start[0] = 0;
1.1326 + count[0] = 1; /* time */
1.1327 + start[1] = hr;
1.1328 + count[1] = 1; /* hour */
1.1329 + start[(ndims-2)] = 0;
1.1330 + count[(ndims-2)] = d1;
1.1331 + start[(ndims-1)] = 0;
1.1332 + count[(ndims-1)] = d2;
1.1333 + if (ndims == 5) {
1.1334 + start[2] = level;
1.1335 + count[2] = 1;
1.1336 + }
1.1337 + switch (nctype) {
1.1338 + case NC_FLOAT:
1.1339 + for (i = 0; i < (d1 * d2); i++)
1.1340 + fvar[i] = (float)var_hours[hr][i];
1.1341 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
1.1342 + break;
1.1343 + case NC_DOUBLE:
1.1344 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
1.1345 + break;
1.1346 + default:
1.1347 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
1.1348 + exit(3);
1.1349 + }
1.1350 + }
1.1351 +
1.1352 + if (nctype == NC_FLOAT)
1.1353 + free(fvar);
1.1354 +
1.1355 + return;
1.1356 +}
1.1357 +
1.1358 +void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
1.1359 +{
1.1360 + int j, k, nlevels;
1.1361 + size_t d1, d2, **cell_samples;
1.1362 + double **means, **stddevs;
1.1363 + struct var *in_vnode, *out_vnode;
1.1364 +
1.1365 + if (input_var_head) {
1.1366 + in_vnode = input_var_head;
1.1367 + /* March through non-metadata variables to compute stats */
1.1368 + for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
1.1369 + if (!in_vnode->metadata) {
1.1370 + /* Make sure time is really the first dimension */
1.1371 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
1.1372 + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
1.1373 + exit(3);
1.1374 + }
1.1375 + /* Figure out input dimensions */
1.1376 + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
1.1377 + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
1.1378 + exit(3);
1.1379 + }
1.1380 + /* Assume 2D output; find dimensions */
1.1381 + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
1.1382 + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
1.1383 + if (in_vnode->ndims == 3)
1.1384 + nlevels = 1;
1.1385 + else
1.1386 + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
1.1387 + /* Allocate working space for means and stddevs */
1.1388 + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
1.1389 + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
1.1390 + printf("Computing statistics for %s\n",
1.1391 + in_vnode->name);
1.1392 + /* Compute means and stddevs, then write them */
1.1393 + for (k = 0; k < nlevels; k++) {
1.1394 + /* Read and compute means */
1.1395 + 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);
1.1396 + /* Find corresponding output variable on the mean output file */
1.1397 + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
1.1398 + /* Write out the means for this variable */
1.1399 + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
1.1400 + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
1.1401 + reset_cell_samples(d1, d2, cell_samples);
1.1402 + /* Read and compute stddevs using means */
1.1403 + 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);
1.1404 + /* Find corresponding output variable on the stddev output file */
1.1405 + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
1.1406 + /* Write out stddevs for this variable */
1.1407 + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
1.1408 + }
1.1409 +
1.1410 + /* Free working space for means and stddevs */
1.1411 + free_means_stddevs(means, stddevs, cell_samples);
1.1412 + }
1.1413 + in_vnode = in_vnode->next;
1.1414 + }
1.1415 + }
1.1416 + return;
1.1417 +}
1.1418 +
1.1419 +void usage(char *program)
1.1420 +{
1.1421 + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
1.1422 + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
1.1423 + exit(4);
1.1424 +}
1.1425 +
1.1426 +int main(int argc, char **argv)
1.1427 +{
1.1428 + int i, ic, nifnames;
1.1429 + size_t len, pos, nsamples;
1.1430 + char *mfname, *sfname, **ifnames, *flist;
1.1431 +
1.1432 + ifnames = NULL;
1.1433 + mfname = sfname = flist = NULL;
1.1434 + nifnames = 0;
1.1435 +
1.1436 +#ifdef DEBUG
1.1437 + printf("Number of metadata variables (nmvars) = %d\n", nmvars);
1.1438 + fflush(stdout);
1.1439 +#endif /* DEBUG */
1.1440 +
1.1441 + /* check command-line flags and switches */
1.1442 + while ((ic = getopt(argc, argv, "m:s:")) != -1) {
1.1443 + switch(ic) {
1.1444 + case 'm':
1.1445 + if (!(mfname = strdup(optarg))) {
1.1446 + perror("mfname");
1.1447 + exit(2);
1.1448 + }
1.1449 + break;
1.1450 + case 's':
1.1451 + if (!(sfname = strdup(optarg))) {
1.1452 + perror("sfname");
1.1453 + exit(2);
1.1454 + }
1.1455 + break;
1.1456 + }
1.1457 + }
1.1458 +
1.1459 + if (!mfname) {
1.1460 + fprintf(stderr, "Output file name for writing means required.\n");
1.1461 + usage(argv[0]);
1.1462 + }
1.1463 + if (!sfname) {
1.1464 + fprintf(stderr, "Output file name for writing standard deviations required.\n");
1.1465 + usage(argv[0]);
1.1466 + }
1.1467 +
1.1468 + if (optind < argc) {
1.1469 + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
1.1470 + perror("ifnames");
1.1471 + exit(2);
1.1472 + }
1.1473 + for (i = optind; i < argc; i++) {
1.1474 + if (!(ifnames[nifnames++] = strdup(argv[i]))) {
1.1475 + perror("ifnames[nifnames++]");
1.1476 + exit(2);
1.1477 + }
1.1478 + }
1.1479 + }
1.1480 + else {
1.1481 + fprintf(stderr, "At least one input file name is required.\n");
1.1482 + usage(argv[0]);
1.1483 + }
1.1484 +
1.1485 +
1.1486 + /*
1.1487 + * Build list of input files to be included in the global history
1.1488 + * attribute on the two outputs files.
1.1489 + */
1.1490 + for (i = len = 0; i < nifnames; i++)
1.1491 + len += strlen(ifnames[i]);
1.1492 + len += nifnames + 1; /* space in front of every name + null terminator */
1.1493 + if (!(flist = (char *)malloc(len * sizeof(char)))) {
1.1494 + perror("flist");
1.1495 + exit(2);
1.1496 + }
1.1497 + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
1.1498 + sprintf(&flist[pos], " %s", ifnames[i]);
1.1499 +#ifdef DEBUG
1.1500 + printf("flist=%s\n", flist);
1.1501 + fflush(stdout);
1.1502 +#endif /* DEBUG */
1.1503 +
1.1504 + /* Open every file to count up the number of time samples. */
1.1505 + nsamples = count_samples(nifnames, ifnames);
1.1506 +#ifdef DEBUG
1.1507 + printf("Number of samples across all files = %ld\n", (long)nsamples);
1.1508 + fflush(stdout);
1.1509 +#endif /* DEBUG */
1.1510 +
1.1511 + /*
1.1512 + * Open the *last* input file and the two output files (for mean and
1.1513 + * standard deviation). Define dimensions, variables, and attributes
1.1514 + * for the two output files based on the *last* input files. The
1.1515 + * last file is used because some metadata variables will contain
1.1516 + * only values from the last time slice from the period over which
1.1517 + * the statistics are computed.
1.1518 + */
1.1519 + open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples);
1.1520 +
1.1521 + compute_stats(nifnames, ifnames, nsamples);
1.1522 +
1.1523 + /* Close the two output files */
1.1524 + wrap_nc(nc_close(mean_ncid));
1.1525 + wrap_nc(nc_close(stddev_ncid));
1.1526 +
1.1527 + return 0;
1.1528 +}