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