1.1 --- a/Makefile Wed Oct 03 11:23:02 2007 -0400
1.2 +++ b/Makefile Wed Oct 10 11:59:02 2007 -0400
1.3 @@ -15,7 +15,7 @@
1.4 # robin1 and Penguins
1.5 CFLAGS=-g -Wall -O $(CPPFLAGS)
1.6
1.7 -all: h1_summary h1_summary2 add_total_fields
1.8 +all: h1_summary h1_summary2 h1_summary_cb add_total_fields
1.9
1.10 h1_summary: h1_summary.o
1.11 $(CC) $(CFLAGS) -o $@ h1_summary.o $(LIBS)
1.12 @@ -23,16 +23,21 @@
1.13 h1_summary2: h1_summary2.o
1.14 $(CC) $(CFLAGS) -o $@ h1_summary2.o $(LIBS)
1.15
1.16 +h1_summary_cb: h1_summary_cb.o
1.17 + $(CC) $(CFLAGS) -o $@ h1_summary_cb.o $(LIBS)
1.18 +
1.19 add_total_fields: add_total_fields.o
1.20 $(CC) $(CFLAGS) -o $@ add_total_fields.o $(LIBS)
1.21
1.22 clean:
1.23 $(RM) -f h1_summary.o h1_summary
1.24 $(RM) -f h1_summary2.o h1_summary2
1.25 + $(RM) -f h1_summary_cb.o h1_summary_cb
1.26 $(RM) -f add_total_fields.o add_total_fields
1.27
1.28 install: all
1.29 cp -p h1_summary $(HOME)/bin/h1_summary.`uname -n`
1.30 cp -p h1_summary2 $(HOME)/bin/h1_summary2.`uname -n`
1.31 + cp -p h1_summary_cb $(HOME)/bin/h1_summary_cb.`uname -n`
1.32 cp -p add_total_fields $(HOME)/bin/add_total_fields.`uname -n`
1.33
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/h1_summary_cb.c Wed Oct 10 11:59:02 2007 -0400
2.3 @@ -0,0 +1,1737 @@
2.4 +#include <stdio.h>
2.5 +#include <unistd.h>
2.6 +#include <stdlib.h>
2.7 +#include <string.h>
2.8 +#include <math.h>
2.9 +#include "netcdf.h"
2.10 +
2.11 +/*
2.12 + * Loop through one month of hourly history tapes from the Community Land Model
2.13 + * and generate monthly statistics (means and standard deviations) of fields
2.14 + * by hour of day.
2.15 + */
2.16 +
2.17 +#define MEAN_TYPE 0
2.18 +#define STDDEV_TYPE 1
2.19 +#define NUM_TYPES 2
2.20 +
2.21 +#define HOURS_PER_DAY 24
2.22 +#define MIN_PER_HOUR 60
2.23 +#define SEC_PER_MIN 60
2.24 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
2.25 +
2.26 +static char *metadata_vars[] = {
2.27 + "lon",
2.28 + "lat",
2.29 + "lonrof",
2.30 + "latrof",
2.31 + "time",
2.32 + "hour", /* new metadata variable to be added to output files */
2.33 + "levsoi",
2.34 + "levlak",
2.35 + "edgen",
2.36 + "edgee",
2.37 + "edges",
2.38 + "edgew",
2.39 + "longxy",
2.40 + "latixy",
2.41 + "area",
2.42 + "areaupsc",
2.43 + "topo",
2.44 + "topodnsc",
2.45 + "landfrac",
2.46 + "landmask",
2.47 + "pftmask",
2.48 + "indxupsc",
2.49 + "mcdate",
2.50 + "mcsec",
2.51 + "mdcur",
2.52 + "mscur",
2.53 + "nstep",
2.54 + "time_bounds",
2.55 + "date_written",
2.56 + "time_written"
2.57 +};
2.58 +
2.59 +struct text_att {
2.60 + char *name;
2.61 + char *value;
2.62 + struct text_att *next;
2.63 + struct text_att *prev;
2.64 +};
2.65 +
2.66 +struct dim {
2.67 + int dimid;
2.68 + char *name;
2.69 + size_t len;
2.70 + struct dim *next;
2.71 + struct dim *prev;
2.72 +};
2.73 +
2.74 +struct var {
2.75 + int ncvarid;
2.76 + char *name;
2.77 + nc_type nctype;
2.78 + int ndims;
2.79 + int *dimids;
2.80 + int natts;
2.81 + char metadata;
2.82 + char FillFlag;
2.83 + float FillValue;
2.84 + struct var *next;
2.85 + struct var *prev;
2.86 +};
2.87 +
2.88 +static char *time_name = "time";
2.89 +static char *time_bounds_attname = "bounds";
2.90 +static char *time_bounds_name = "time_bounds";
2.91 +static char *climatology_bounds_attname = "climatology";
2.92 +static char *climatology_bounds_name = "climatology_bounds";
2.93 +static char *mcdate_name = "mcdate";
2.94 +static char *mcsec_name = "mcsec";
2.95 +static char *history_name = "history";
2.96 +static char *nsamples_name = "nsamples";
2.97 +static char *cell_method_name = "cell_method";
2.98 +static char *cell_methods_prefix[] = { "time: mean within days time: mean over days", "time: mean within days times: standard_deviation over days" };
2.99 +/* input stuff */
2.100 +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
2.101 +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
2.102 +static struct text_att *input_gatt_head = NULL;
2.103 +static struct dim *input_dim_head = NULL, **input_dim_idx;
2.104 +static struct var *input_var_head = NULL;
2.105 +/* translation stuff */
2.106 +static int *idid2mdid, *idid2sdid; /* dimension IDs */
2.107 +/* output stuff */
2.108 +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
2.109 +static struct dim *mean_dim_head = NULL;
2.110 +static struct var *mean_var_head = NULL;
2.111 +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
2.112 +static struct dim *stddev_dim_head = NULL;
2.113 +static struct var *stddev_var_head = NULL;
2.114 +/* output values */
2.115 +static double *times;
2.116 +static double *tbounds;
2.117 +static int *mcdate;
2.118 +static int *mcsec;
2.119 +
2.120 +char is_metadata(char *name)
2.121 +{
2.122 + int i;
2.123 +
2.124 + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
2.125 +
2.126 + if (i < nmvars)
2.127 + return 1;
2.128 + else
2.129 + return 0;
2.130 +}
2.131 +
2.132 +#if 0
2.133 +struct dim *dimlist_find_by_name(struct dim *head, char *name)
2.134 +{
2.135 + int i;
2.136 + struct dim *p = NULL;
2.137 +
2.138 + if (head) {
2.139 + node = head;
2.140 + for (i = 0 ; i == 0 || node != head; i++) {
2.141 + if (!strcmp(node->name, name))
2.142 + p = node;
2.143 + node = node->next;
2.144 + }
2.145 + return p;
2.146 + }
2.147 + else
2.148 + return NULL;
2.149 +}
2.150 +#endif
2.151 +
2.152 +struct var *varlist_find_by_name(struct var *head, char *name)
2.153 +{
2.154 + int i;
2.155 + struct var *node;
2.156 +
2.157 + if (head) {
2.158 + node = head;
2.159 + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
2.160 + node = node->next;
2.161 + if (i && node == head)
2.162 + return NULL;
2.163 + else
2.164 + return node;
2.165 + }
2.166 + else
2.167 + return NULL;
2.168 +}
2.169 +
2.170 +void gattlist_add_node(struct text_att **headp, char *name, char *value)
2.171 +{
2.172 + struct text_att *head, *node;
2.173 +
2.174 + head = *headp;
2.175 +
2.176 + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
2.177 + perror("gattlist_add_node");
2.178 + exit(2);
2.179 + }
2.180 +
2.181 + if (!(node->name = strdup(name))) {
2.182 + perror("gattlist_add_node");
2.183 + exit(2);
2.184 + }
2.185 + if (!(node->value = strdup(value))) {
2.186 + perror("gattlist_add_node");
2.187 + exit(2);
2.188 + }
2.189 +
2.190 + if (head == NULL) {
2.191 + /* first one */
2.192 + *headp = node;
2.193 + node->next = node;
2.194 + node->prev = node;
2.195 + }
2.196 + else {
2.197 + node->next = head;
2.198 + node->prev = head->prev;
2.199 + /* set this after setting node->prev from here */
2.200 + head->prev = node;
2.201 + /* set this after having set node->prev */
2.202 + node->prev->next = node;
2.203 + }
2.204 +
2.205 + return;
2.206 +}
2.207 +
2.208 +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
2.209 +{
2.210 + struct dim *head, *node;
2.211 +
2.212 + head = *headp;
2.213 +
2.214 + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
2.215 + perror("dimlist_add_node");
2.216 + exit(2);
2.217 + }
2.218 +
2.219 + node->dimid = dimid;
2.220 + if (!(node->name = strdup(name))) {
2.221 + perror("dimlist_add_node");
2.222 + exit(2);
2.223 + }
2.224 + node->len = len;
2.225 +
2.226 + if (head == NULL) {
2.227 + /* first one */
2.228 + *headp = node;
2.229 + node->next = node;
2.230 + node->prev = node;
2.231 + }
2.232 + else {
2.233 + node->next = head;
2.234 + node->prev = head->prev;
2.235 + /* set this after setting node->prev from here */
2.236 + head->prev = node;
2.237 + /* set this after having set node->prev */
2.238 + node->prev->next = node;
2.239 + }
2.240 +
2.241 + return node;
2.242 +}
2.243 +
2.244 +void varlist_add_node(struct var **headp, int ncvarid, char *name,
2.245 + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
2.246 + float FillValue)
2.247 +{
2.248 + int i;
2.249 + struct var *head, *node;
2.250 +
2.251 + head = *headp;
2.252 +
2.253 + if (!(node = (struct var *)malloc(sizeof(struct var)))) {
2.254 + perror("varlist_add_node");
2.255 + exit(3);
2.256 + }
2.257 +
2.258 + node->ncvarid = ncvarid;
2.259 + if (!(node->name = strdup(name))) {
2.260 + perror("varlist_add_node");
2.261 + exit(3);
2.262 + }
2.263 + node->nctype = nctype;
2.264 + node->ndims = ndims;
2.265 + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
2.266 + perror("varlist_add_node");
2.267 + exit(3);
2.268 + }
2.269 + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
2.270 + node->natts = natts;
2.271 + node->metadata = is_metadata(name);
2.272 + node->FillFlag = FillFlag;
2.273 + node->FillValue = FillValue;
2.274 +
2.275 + if (head == NULL) {
2.276 + /* first one */
2.277 + *headp = node;
2.278 + node->next = node;
2.279 + node->prev = node;
2.280 + }
2.281 + else {
2.282 + node->next = head;
2.283 + node->prev = head->prev;
2.284 + /* set this after setting node->prev from here */
2.285 + head->prev = node;
2.286 + /* set this after having set node->prev */
2.287 + node->prev->next = node;
2.288 + }
2.289 +
2.290 + return;
2.291 +}
2.292 +
2.293 +void varlist_print(struct var *headp)
2.294 +{
2.295 + int i, j;
2.296 + struct var *node;
2.297 +
2.298 + printf("Beginning of Variable List\n");
2.299 +
2.300 + if (headp) {
2.301 + node = headp;
2.302 + for (j = 0; j == 0 || node != headp; j++) {
2.303 + printf("Variable %d (ptr=%ld):\n", j, (long)node);
2.304 + printf("\tncvarid = %d\n", node->ncvarid);
2.305 + printf("\tname = %s\n", node->name);
2.306 + printf("\tnctype = %d\n", node->nctype);
2.307 + printf("\tndims = %d\n", node->ndims);
2.308 + printf("\tdimids =");
2.309 + for (i = 0; i < node->ndims; i++)
2.310 + printf(" %d", node->dimids[i]);
2.311 + printf("\n");
2.312 + printf("\tnatts = %d\n", node->natts);
2.313 + printf("\tmetadata = %d\n", (int)node->metadata);
2.314 + printf("\tnext = %ld\n", (long)node->next);
2.315 + printf("\tprev = %ld\n", (long)node->prev);
2.316 + node = node->next;
2.317 + }
2.318 + }
2.319 + else {
2.320 + printf("\tList undefined\n");
2.321 + }
2.322 +
2.323 + printf("End of Variable List\n");
2.324 +
2.325 + return;
2.326 +}
2.327 +
2.328 +void wrap_nc(int status)
2.329 +{
2.330 + if (status != NC_NOERR) {
2.331 + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
2.332 + exit(1);
2.333 + }
2.334 +
2.335 + return;
2.336 +}
2.337 +
2.338 +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
2.339 +{
2.340 + int i;
2.341 + char name[NC_MAX_NAME+1];
2.342 + size_t len;
2.343 + struct dim **in_dim_idx;
2.344 +
2.345 + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
2.346 + perror("*in_dim_idxp");
2.347 + exit(3);
2.348 + }
2.349 + in_dim_idx = *in_dim_idxp;
2.350 +
2.351 + for (i = 0; i < ndims; i++) {
2.352 + wrap_nc(nc_inq_dim(ncid, i, name, &len));
2.353 + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
2.354 + }
2.355 +
2.356 + return;
2.357 +}
2.358 +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
2.359 +{
2.360 + int i;
2.361 + char name[NC_MAX_NAME+1], *value;
2.362 + nc_type xtype;
2.363 + size_t len;
2.364 +
2.365 + for (i = 0; i < ngatts; i++) {
2.366 + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
2.367 + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
2.368 + if (xtype != NC_CHAR) {
2.369 + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
2.370 + exit(2);
2.371 + }
2.372 + if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
2.373 + perror("get_gatts: value");
2.374 + exit(3);
2.375 + }
2.376 + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
2.377 + value[(len+1-1)] = '\0';
2.378 + gattlist_add_node(gatt_headp, name, value);
2.379 + free(value); /* because gattlist_add_node() duplicates it */
2.380 + }
2.381 +
2.382 + return;
2.383 +}
2.384 +
2.385 +void get_vars(int ncid, int nvars, struct var **var_headp)
2.386 +{
2.387 + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
2.388 + size_t len;
2.389 + float FillValue;
2.390 + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
2.391 + nc_type xtype, atype;
2.392 +
2.393 + for (i = 0; i < nvars; i++) {
2.394 + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
2.395 + &natts));
2.396 + /* Check for _FillValue */
2.397 + FillFlag = 0;
2.398 + FillValue = 0.;
2.399 + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
2.400 + == NC_NOERR) {
2.401 + if (atype == NC_FLOAT && len) {
2.402 + wrap_nc(nc_get_att_float(ncid, i,
2.403 + fill_att_name, &FillValue));
2.404 + FillFlag = 1;
2.405 + }
2.406 + }
2.407 + varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
2.408 + natts, FillFlag, FillValue);
2.409 + }
2.410 +
2.411 + return;
2.412 +}
2.413 +
2.414 +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
2.415 + size_t nsamples, int **in2outp, int out_ncid,
2.416 + struct dim **out_dim_headp, int *out_unlimdimidp)
2.417 +{
2.418 + /*
2.419 + * Define dimensions on new files and build translation tables between
2.420 + * dimension IDs.
2.421 + */
2.422 + int j, dimid, ndims, *in2out;
2.423 + size_t len;
2.424 + struct dim *dnode;
2.425 +
2.426 + ndims = 0;
2.427 +
2.428 + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
2.429 + perror("put_dimensions:in2outp");
2.430 + exit(3);
2.431 + }
2.432 + in2out = *in2outp;
2.433 +
2.434 + if (in_dim_head) {
2.435 + dnode = in_dim_head;
2.436 + for (j = 0; j == 0 || dnode != in_dim_head; j++) {
2.437 + if (dnode->dimid == in_unlimdimid)
2.438 + len = NC_UNLIMITED;
2.439 + else
2.440 + len = dnode->len;
2.441 + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
2.442 + dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
2.443 + ++ndims;
2.444 + in2out[dnode->dimid] = dimid;
2.445 + if (dnode->dimid == in_unlimdimid)
2.446 + *out_unlimdimidp = dimid;
2.447 + /* Just after the "time" dimension, add the new
2.448 + * "nsamples" dimension. */
2.449 + if (!strcmp(dnode->name, time_name)) {
2.450 + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
2.451 + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
2.452 + ++ndims;
2.453 + }
2.454 + dnode = dnode->next;
2.455 + }
2.456 + }
2.457 + else {
2.458 + fprintf(stderr, "WARNING: No dimensions defined!\n");
2.459 + }
2.460 +
2.461 + return ndims;
2.462 +}
2.463 +
2.464 +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
2.465 +{
2.466 + /*
2.467 + * Write out global attributes matching those of the input file.
2.468 + * Change history attribute to the string passed in as an argument.
2.469 + */
2.470 + int j, hflag = 0, ngatts;
2.471 + char *value;
2.472 + struct text_att *anode;
2.473 +
2.474 + ngatts = 0;
2.475 +
2.476 + if (in_gatt_head) {
2.477 + anode = in_gatt_head;
2.478 + for (j = 0; j == 0 || anode != in_gatt_head; j++) {
2.479 + if (!strcmp(anode->name, history_name)) {
2.480 + value = out_history;
2.481 + ++hflag;
2.482 + }
2.483 + else
2.484 + value = anode->value;
2.485 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
2.486 + ++ngatts;
2.487 + anode = anode->next;
2.488 + }
2.489 + /* If no history attribute on input, add one to the output */
2.490 + if (!hflag) {
2.491 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
2.492 + ++ngatts;
2.493 + }
2.494 + }
2.495 + else {
2.496 + fprintf(stderr, "WARNING: No global attributes defined!\n");
2.497 + }
2.498 +
2.499 + return ngatts;
2.500 +}
2.501 +
2.502 +int translate_dimids(int in_ndims, int *in_dimids, int *in2out, int *out_dimids)
2.503 +{
2.504 + /*
2.505 + * Translate between input dimension IDs and those for the output file.
2.506 + */
2.507 + int i, ndims;
2.508 +
2.509 + for (i = ndims = 0; i < in_ndims; i++, ndims++)
2.510 + out_dimids[ndims] = in2out[in_dimids[i]];
2.511 +
2.512 + return ndims;
2.513 +}
2.514 +
2.515 +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
2.516 + int in_varid, int out_ncid, int out_varid, char **cell_methods)
2.517 +{
2.518 + /*
2.519 + * Copy the attributes of the input variable to those of the output
2.520 + * variable. If the variable is not a metadata variable, ensure that
2.521 + * the cell_method attribute is set correctly, depending on the output
2.522 + * file type.
2.523 + */
2.524 +
2.525 + int i, natts;
2.526 + char name[NC_MAX_NAME+1], cmflag = 0;
2.527 +
2.528 + for (i = natts = 0; i < in_natts; i++, natts++) {
2.529 + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
2.530 + if (!strcmp(name, cell_method_name) && !metadata) {
2.531 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
2.532 + cmflag = 1;
2.533 + }
2.534 + else if (!strcmp(name, time_bounds_attname))
2.535 + /* Change time_bounds to climatology_bounds */
2.536 + wrap_nc(nc_put_att_text(out_ncid, out_varid, climatology_bounds_attname, strlen(climatology_bounds_name), climatology_bounds_name));
2.537 + else
2.538 + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
2.539 + }
2.540 + /*
2.541 + * If no cell_method attribute was specified for a non-metadata
2.542 + * variable on the input file, add the appropriate cell_method anyway
2.543 + * on the output file.
2.544 + */
2.545 + if (!cmflag && !metadata) {
2.546 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
2.547 + natts++;
2.548 + }
2.549 +
2.550 + return natts;
2.551 +}
2.552 +
2.553 +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
2.554 + int *in2out, int out_ncid, struct var **out_var_headp,
2.555 + char **cell_methods, int stat_type)
2.556 +{
2.557 + /*
2.558 + * Define variables on output file and copy attributes from input file.
2.559 + * Return number of variables defined.
2.560 + */
2.561 + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
2.562 + struct var *vnode;
2.563 +
2.564 + nvars = 0;
2.565 +
2.566 + if (in_var_head) {
2.567 + vnode = in_var_head;
2.568 + /*
2.569 + * March through input variables creating (mostly) the same
2.570 + * variables on the output file.
2.571 + */
2.572 + for (j = 0; j == 0 || vnode != in_var_head; j++) {
2.573 + ndims = translate_dimids(vnode->ndims, vnode->dimids, in2out, dimids);
2.574 + /* Create all normal data variables, but only the four
2.575 + * time-based metadata variables: time,
2.576 + * climatology_bounds, mcdate, and mcsec. */
2.577 + if (!strcmp(vnode->name, time_name)) {
2.578 + /* Magically promote "time" variable to
2.579 + * double precision on output file */
2.580 + wrap_nc(nc_def_var(out_ncid, vnode->name, NC_DOUBLE, ndims, dimids, &varid));
2.581 + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);
2.582 + varlist_add_node(out_var_headp, varid, vnode->name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
2.583 + ++nvars;
2.584 + /* Create "climatology_bounds" right after
2.585 + * time and instead of "time_bounds", with no
2.586 + * attributes */
2.587 + wrap_nc(nc_def_var(out_ncid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, &varid));
2.588 + natts = 0;
2.589 + varlist_add_node(out_var_headp, varid, climatology_bounds_name, NC_DOUBLE, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
2.590 + ++nvars;
2.591 + }
2.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))) {
2.593 + /* If it is not a metadata variable OR it is
2.594 + * mcdate OR it is mcsec OR (it is a metadata
2.595 + * variable that does not have time as its
2.596 + * first dimension), then create it */
2.597 + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
2.598 + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid, cell_methods);
2.599 + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
2.600 + ++nvars;
2.601 + }
2.602 + vnode = vnode->next;
2.603 + }
2.604 + }
2.605 + else {
2.606 + fprintf(stderr, "WARNING: No variables defined!\n");
2.607 + }
2.608 +
2.609 + return nvars;
2.610 +}
2.611 +
2.612 +void *alloc_var(nc_type nctype, size_t len)
2.613 +{
2.614 + void *val;
2.615 +
2.616 + switch (nctype) {
2.617 + case NC_FLOAT:
2.618 + if (!(val = (float *)malloc((len) * sizeof(float)))) {
2.619 + perror("alloc_var: val");
2.620 + exit(3);
2.621 + }
2.622 + break;
2.623 + case NC_INT:
2.624 + if (!(val = (int *)malloc((len) * sizeof(int)))) {
2.625 + perror("alloc_var: val");
2.626 + exit(3);
2.627 + }
2.628 + break;
2.629 + case NC_DOUBLE:
2.630 + if (!(val = (double *)malloc((len) * sizeof(double)))) {
2.631 + perror("alloc_var: val");
2.632 + exit(3);
2.633 + }
2.634 + break;
2.635 + case NC_CHAR:
2.636 + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
2.637 + perror("alloc_var: val");
2.638 + exit(3);
2.639 + }
2.640 + break;
2.641 + default:
2.642 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.643 + exit(3);
2.644 + }
2.645 +
2.646 + return val;
2.647 +}
2.648 +
2.649 +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
2.650 +{
2.651 + int i;
2.652 + size_t len = (size_t)1;
2.653 + void *val;
2.654 +
2.655 + /* Figure out the total size */
2.656 + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
2.657 +
2.658 + val = alloc_var(nctype, len);
2.659 +
2.660 + switch (nctype) {
2.661 + case NC_FLOAT:
2.662 + wrap_nc(nc_get_var_float(ncid, varid, val));
2.663 + break;
2.664 + case NC_INT:
2.665 + wrap_nc(nc_get_var_int(ncid, varid, val));
2.666 + break;
2.667 + case NC_DOUBLE:
2.668 + wrap_nc(nc_get_var_double(ncid, varid, val));
2.669 + break;
2.670 + case NC_CHAR:
2.671 + wrap_nc(nc_get_var_text(ncid, varid, val));
2.672 + break;
2.673 + default:
2.674 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.675 + exit(3);
2.676 + }
2.677 +
2.678 + return val;
2.679 +}
2.680 +
2.681 +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
2.682 +{
2.683 + int i;
2.684 + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.685 + void *val;
2.686 +
2.687 + /* Make sure time is really the first dimension */
2.688 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
2.689 + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
2.690 + exit(3);
2.691 + }
2.692 + /*
2.693 + * Figure out the total size for the slice (start at second dimension)
2.694 + * and build start/count vectors.
2.695 + */
2.696 + start[0] = tslice;
2.697 + count[0] = 1;
2.698 + for (i = 1; i < ndims; i++) {
2.699 + start[i] = 0;
2.700 + count[i] = (dim_idx[dimids[i]])->len;
2.701 + len *= (dim_idx[dimids[i]])->len;
2.702 + }
2.703 +
2.704 + val = alloc_var(nctype, len);
2.705 +
2.706 + switch (nctype) {
2.707 + case NC_FLOAT:
2.708 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
2.709 + break;
2.710 + case NC_INT:
2.711 + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
2.712 + break;
2.713 + case NC_DOUBLE:
2.714 + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
2.715 + break;
2.716 + case NC_CHAR:
2.717 + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
2.718 + break;
2.719 + default:
2.720 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.721 + exit(3);
2.722 + }
2.723 +
2.724 + return val;
2.725 +}
2.726 +
2.727 +void write_var(int ncid, int varid, nc_type nctype, void *val)
2.728 +{
2.729 + switch (nctype) {
2.730 + case NC_FLOAT:
2.731 + wrap_nc(nc_put_var_float(ncid, varid, val));
2.732 + break;
2.733 + case NC_INT:
2.734 + wrap_nc(nc_put_var_int(ncid, varid, val));
2.735 + break;
2.736 + case NC_DOUBLE:
2.737 + wrap_nc(nc_put_var_double(ncid, varid, val));
2.738 + break;
2.739 + case NC_CHAR:
2.740 + wrap_nc(nc_put_var_text(ncid, varid, val));
2.741 + break;
2.742 + default:
2.743 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.744 + exit(3);
2.745 + }
2.746 +
2.747 + return;
2.748 +}
2.749 +
2.750 +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
2.751 +{
2.752 + int i;
2.753 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.754 +
2.755 + /* Make sure time is really the first dimension */
2.756 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
2.757 + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
2.758 + exit(3);
2.759 + }
2.760 +
2.761 + /* Build start/count vectors */
2.762 + start[0] = tslice;
2.763 + count[0] = 1;
2.764 + for (i = 1; i < ndims; i++) {
2.765 + start[i] = 0;
2.766 + count[i] = (dim_idx[dimids[i]])->len;
2.767 + }
2.768 +
2.769 + switch (nctype) {
2.770 + case NC_FLOAT:
2.771 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
2.772 + break;
2.773 + case NC_INT:
2.774 + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
2.775 + break;
2.776 + case NC_DOUBLE:
2.777 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
2.778 + break;
2.779 + case NC_CHAR:
2.780 + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
2.781 + break;
2.782 + default:
2.783 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.784 + exit(3);
2.785 + }
2.786 +
2.787 + return;
2.788 +}
2.789 +
2.790 +void write_time_var(int ncid, int varid, nc_type nctype, void *time_var)
2.791 +{
2.792 + size_t start[1], count[1];
2.793 +
2.794 + start[0] = 0;
2.795 + count[0] = 24;
2.796 +
2.797 + switch (nctype) {
2.798 + case NC_FLOAT:
2.799 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, time_var));
2.800 + break;
2.801 + case NC_INT:
2.802 + wrap_nc(nc_put_vara_int(ncid, varid, start, count, time_var));
2.803 + break;
2.804 + case NC_DOUBLE:
2.805 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, time_var));
2.806 + break;
2.807 + case NC_CHAR:
2.808 + wrap_nc(nc_put_vara_text(ncid, varid, start, count, time_var));
2.809 + break;
2.810 + default:
2.811 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.812 + exit(3);
2.813 + }
2.814 +
2.815 + return;
2.816 +}
2.817 +
2.818 +void write_tbounds(int ncid, int varid, double *tbounds)
2.819 +{
2.820 + size_t start[2], count[2];
2.821 +
2.822 + start[0] = start[1] = 0;
2.823 + count[0] = 24; count[1] = 2;
2.824 +
2.825 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, tbounds));
2.826 +
2.827 + return;
2.828 +}
2.829 +
2.830 +void copy_metadata(int in_ncid, struct var *in_var_head,
2.831 + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head,
2.832 + double *times, double *tbounds, int *mcdate, int *mcsec)
2.833 +{
2.834 + int i;
2.835 + struct var *in_vnode, *out_vnode;
2.836 + void *val;
2.837 +
2.838 + if (!in_var_head) return;
2.839 +
2.840 + in_vnode = in_var_head;
2.841 + /*
2.842 + * March through input variables to stuff metadata values into
2.843 + * the new output files. NOTE: The only time-based metadata
2.844 + * variables, explicitly handled, are: time, climatology_bounds,
2.845 + * mcdate, and mcsec. These are passed in.
2.846 + */
2.847 + for (i = 0; i == 0 || in_vnode != in_var_head; i++) {
2.848 + if (in_vnode->metadata) {
2.849 + /* Find the corresponding output variable */
2.850 + if (!strcmp(in_vnode->name, time_bounds_name))
2.851 + out_vnode = varlist_find_by_name(out_var_head, climatology_bounds_name);
2.852 + else
2.853 + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
2.854 + if (out_vnode) {
2.855 + /* Read values from input and write to output */
2.856 + if (strcmp((in_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
2.857 + /* "time" is not the first dimension,
2.858 + * so just copy the variable values */
2.859 +#ifdef DEBUG
2.860 + printf("Copying metadata variable %s\n",
2.861 + in_vnode->name);
2.862 +#endif /* DEBUG */
2.863 + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
2.864 + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
2.865 + free(val);
2.866 + }
2.867 + else {
2.868 + /* "time" is the first dimension, so
2.869 + * explicitly handle the "time",
2.870 + * "climatology_bounds", "mcdate", and
2.871 + * "mcsec" while dropping all others */
2.872 + if (!strcmp(out_vnode->name, time_name)) {
2.873 + /* Force time stamps to be
2.874 + * those computed previously,
2.875 + * and remember these are now
2.876 + * double precision. */
2.877 +#ifdef DEBUG
2.878 + printf("Writing metadata variable %s\n", out_vnode->name);
2.879 +#endif /* DEBUG */
2.880 + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, times);
2.881 + }
2.882 + if (!strcmp(out_vnode->name, climatology_bounds_name)) {
2.883 + /* time_bounds are now
2.884 + * climatology_bounds and have
2.885 + * pre-computed values */
2.886 +#ifdef DEBUG
2.887 + printf("Writing metadata variable %s\n", out_vnode->name);
2.888 +#endif /* DEBUG */
2.889 + write_tbounds(out_ncid, out_vnode->ncvarid, tbounds);
2.890 + }
2.891 + else if (!strcmp(out_vnode->name, mcdate_name)) {
2.892 + /* mcdate is pre-computed */
2.893 +#ifdef DEBUG
2.894 + printf("Writing metadata variable %s\n", out_vnode->name);
2.895 +#endif /* DEBUG */
2.896 + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcdate);
2.897 + }
2.898 + else if (!strcmp(out_vnode->name, mcsec_name)) {
2.899 + /* mcsec is pre-computed */
2.900 +#ifdef DEBUG
2.901 + printf("Writing metadata variable %s\n", out_vnode->name);
2.902 +#endif /* DEBUG */
2.903 + write_time_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, mcsec);
2.904 + }
2.905 + }
2.906 + }
2.907 + else {
2.908 +#ifdef DEBUG
2.909 + printf("Dropping metadata variable %s\n", in_vnode->name);
2.910 +#endif /* DEBUG */
2.911 + }
2.912 + }
2.913 + else {
2.914 +#ifdef DEBUG
2.915 + printf("Skipping non-metadata variable %s\n",
2.916 + in_vnode->name);
2.917 +#endif /* DEBUG */
2.918 + }
2.919 +#ifdef DEBUG
2.920 + printf ("Done\n");
2.921 +#endif /* DEBUG */
2.922 + in_vnode = in_vnode->next;
2.923 + }
2.924 +
2.925 + return;
2.926 +}
2.927 +
2.928 +void get_time_bounds(char *first_fname, char *last_fname, double *times, double *tbounds, int *mcdate, int *mcsec)
2.929 +{
2.930 + int i, time_dimid, time_bounds_varid, mcdate_varid, mcsec_varid;
2.931 + size_t time_len, time_bounds_start[2], time_bounds_count[2],
2.932 + mc_start[1], mc_count[1];
2.933 + double bnd1[24], bnd2[24], day1, day2;
2.934 +
2.935 + /* Open first input file */
2.936 + wrap_nc(nc_open(first_fname, NC_NOWRITE, &input_ncid));
2.937 + /* Get dimension ID of the time dimension */
2.938 + wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
2.939 + /* Get length of time dimension */
2.940 + wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
2.941 + /* Get variable ID of time_bounds variable */
2.942 + wrap_nc(nc_inq_varid(input_ncid, time_bounds_name, &time_bounds_varid));
2.943 + /* Set start/count for reading the left values of time_bounds from the
2.944 + * first 24 timeslices of the first file */
2.945 + time_bounds_start[0] = time_bounds_start[1] = 0;
2.946 + time_bounds_count[0] = 24; time_bounds_count[1] = 1;
2.947 + /* Read the values */
2.948 + wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,
2.949 + time_bounds_start, time_bounds_count, bnd1));
2.950 + /* If the first and last file are not the same, close the first one and
2.951 + * open the second one */
2.952 + if (strcmp(first_fname, last_fname)) {
2.953 + /* Close the first input file */
2.954 + wrap_nc(nc_close(input_ncid));
2.955 + /* Open the last input file */
2.956 + wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
2.957 + /* Get dimension ID of the time dimension */
2.958 + wrap_nc(nc_inq_dimid(input_ncid, time_name, &time_dimid));
2.959 + /* Get length of time dimension */
2.960 + wrap_nc(nc_inq_dimlen(input_ncid, time_dimid, &time_len));
2.961 + /* Get variable ID of time_bounds variable */
2.962 + wrap_nc(nc_inq_varid(input_ncid, time_bounds_name,
2.963 + &time_bounds_varid));
2.964 + }
2.965 + /* Set start/count for reading the right values of time_bounds from the
2.966 + * last 24 timeslices of the last file */
2.967 + time_bounds_start[0] = time_len - 24; time_bounds_start[1] = 1;
2.968 + time_bounds_count[0] = 24; time_bounds_count[1] = 1;
2.969 + /* Read the value */
2.970 + wrap_nc(nc_get_vara_double(input_ncid, time_bounds_varid,
2.971 + time_bounds_start, time_bounds_count, bnd2));
2.972 + /* Read the last 24 mcdate and mcsec values */
2.973 + wrap_nc(nc_inq_varid(input_ncid, mcdate_name, &mcdate_varid));
2.974 + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &mcsec_varid));
2.975 + mc_start[0] = time_len - 24;
2.976 + mc_count[0] = 24;
2.977 + wrap_nc(nc_get_vara_int(input_ncid, mcdate_varid, mc_start, mc_count,
2.978 + mcdate));
2.979 + wrap_nc(nc_get_vara_int(input_ncid, mcsec_varid, mc_start, mc_count,
2.980 + mcsec));
2.981 +
2.982 + /* Close the last input file */
2.983 + wrap_nc(nc_close(input_ncid));
2.984 +
2.985 + /* Store time bounds and compute new time stamps as midpoints */
2.986 + for (i = 0; i < 24; i++) {
2.987 + tbounds[(i*2)] = bnd1[i];
2.988 + tbounds[(i*2+1)] = bnd2[i];
2.989 + /* Pick one (integer) day near the midpoint of the bounds and
2.990 + * add in the fraction from the right bounds, so for January
2.991 + * the first hour-of-day value would be 15.0 (plus some year
2.992 + * offset) */
2.993 + day1 = floor(bnd1[i]);
2.994 + day2 = floor(bnd2[i]);
2.995 + times[i] = floor((day1 + day2) / 2.0) + (bnd2[i]-day2);
2.996 + }
2.997 +
2.998 + return;
2.999 +}
2.1000 +
2.1001 +char **make_cell_methods(size_t nsamples)
2.1002 +{
2.1003 + int i;
2.1004 + size_t len;
2.1005 + char **cell_methods;
2.1006 +
2.1007 + if (!(cell_methods = (char **)malloc((NUM_TYPES * sizeof(char *))))) {
2.1008 + perror("make_cell_methods:cell_methods");
2.1009 + exit(2);
2.1010 + }
2.1011 + for (i = 0; i < NUM_TYPES; i++) {
2.1012 + len = strlen(cell_methods_prefix[i]) + 48;
2.1013 + if (!(cell_methods[i] = (char *)malloc((len * sizeof(char))))) {
2.1014 + perror("make_cell_methods:cell_methods[i]");
2.1015 + exit(2);
2.1016 + }
2.1017 + sprintf(cell_methods[i], "%s (comment: %ld samples)",
2.1018 + cell_methods_prefix[i], (long)nsamples);
2.1019 + }
2.1020 + return cell_methods;
2.1021 +}
2.1022 +
2.1023 +void open_inout(char *first_fname, char *last_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
2.1024 +{
2.1025 + int i;
2.1026 + char *mean_history_gatt, *stddev_history_gatt,
2.1027 + *mean_prefix="Statistical means from history files:",
2.1028 + *stddev_prefix="Statistical standard deviations from history files:",
2.1029 + **cell_methods;
2.1030 +
2.1031 + /*
2.1032 + * Construct strings for history global attributes for the two output
2.1033 + * files.
2.1034 + */
2.1035 + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
2.1036 + perror("open_inout:mean_history_gatt");
2.1037 + exit(2);
2.1038 + }
2.1039 + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
2.1040 + perror("open_inout:stddev_history_gatt");
2.1041 + exit(2);
2.1042 + }
2.1043 + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
2.1044 + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
2.1045 +
2.1046 + /* Allocate space for time values for each hour of the day */
2.1047 + if (!(times = (double *)malloc((24 * sizeof(double))))) {
2.1048 + perror("open_inout:times");
2.1049 + exit(2);
2.1050 + }
2.1051 + /* Allocate space for two time bounds values for each hour of the day */
2.1052 + if (!(tbounds = (double *)malloc((24 * 2 * sizeof(double))))) {
2.1053 + perror("open_inout:tbounds");
2.1054 + exit(2);
2.1055 + }
2.1056 + /* Allocate space for current date and seconds for each hour of day */
2.1057 + if (!(mcdate = (int *)malloc((24 * sizeof(int))))) {
2.1058 + perror("open_inout:mcdate");
2.1059 + exit(2);
2.1060 + }
2.1061 + if (!(mcsec = (int *)malloc((24 * sizeof(int))))) {
2.1062 + perror("open_inout:mcsec");
2.1063 + exit(2);
2.1064 + }
2.1065 +
2.1066 +
2.1067 + /*
2.1068 + * Explicitly handle "time", "time_bounds" -> "climatology_bounds",
2.1069 + * mcdate, and mcsec metadata variables by constructing their values
2.1070 + * before doing anything else.
2.1071 + */
2.1072 + get_time_bounds(first_fname, last_fname, times, tbounds, mcdate, mcsec);
2.1073 +#ifdef DEBUG
2.1074 + for (i = 0; i < 24; i++)
2.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]);
2.1076 +#endif /* DEBUG */
2.1077 +
2.1078 + /* Open last input file */
2.1079 + wrap_nc(nc_open(last_fname, NC_NOWRITE, &input_ncid));
2.1080 + /* Inquire about number of dimensions, variables, global attributes
2.1081 + * and the ID of the unlimited dimension
2.1082 + */
2.1083 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
2.1084 + &input_unlimdimid));
2.1085 + /* Grab dimension IDs and lengths from input file */
2.1086 + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
2.1087 + /* Grab desired global attributes from input file */
2.1088 + get_gatts(input_ncid, input_ngatts, &input_gatt_head);
2.1089 + /* Grab variable IDs and attributes from input file */
2.1090 + get_vars(input_ncid, input_nvars, &input_var_head);
2.1091 +#ifdef DEBUG
2.1092 + varlist_print(input_var_head);
2.1093 +#endif /* DEBUG */
2.1094 + /* Create netCDF files for monthly means and standard deviations */
2.1095 + /* Create new netCDF files */
2.1096 + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
2.1097 + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
2.1098 + /* Define dimensions */
2.1099 + mean_ndims = put_dimensions(input_dim_head, input_ndims,
2.1100 + input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
2.1101 + &mean_dim_head, &mean_unlimdimid);
2.1102 + stddev_ndims = put_dimensions(input_dim_head, input_ndims,
2.1103 + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
2.1104 + &stddev_dim_head, &stddev_unlimdimid);
2.1105 + /* Define cell_methods[], including number of samples */
2.1106 + cell_methods = make_cell_methods(nsamples);
2.1107 + /* Define variables and their attributes */
2.1108 + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
2.1109 + idid2mdid, mean_ncid, &mean_var_head, cell_methods,
2.1110 + MEAN_TYPE);
2.1111 + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
2.1112 + idid2sdid, stddev_ncid, &stddev_var_head, cell_methods,
2.1113 + STDDEV_TYPE);
2.1114 +#ifdef DEBUG
2.1115 + printf("Mean variables:\n");
2.1116 + varlist_print(mean_var_head);
2.1117 + printf("Stddev variables:\n");
2.1118 + varlist_print(stddev_var_head);
2.1119 +#endif /* DEBUG */
2.1120 + /* Store global attributes */
2.1121 + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
2.1122 + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
2.1123 + stddev_history_gatt);
2.1124 + /* End define mode */
2.1125 + wrap_nc(nc_enddef(mean_ncid));
2.1126 + wrap_nc(nc_enddef(stddev_ncid));
2.1127 + /* Write out metdata variables that do not depend on "time" */
2.1128 + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
2.1129 + mean_var_head, times, tbounds, mcdate, mcsec);
2.1130 + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
2.1131 + stddev_var_head, times, tbounds, mcdate, mcsec);
2.1132 +
2.1133 + wrap_nc(nc_close(input_ncid));
2.1134 +
2.1135 + return;
2.1136 +}
2.1137 +
2.1138 +size_t count_samples(int nifnames, char **input_fnames)
2.1139 +{
2.1140 + int i;
2.1141 + char name[NC_MAX_NAME+1];
2.1142 + size_t len, cnt = 0;
2.1143 +
2.1144 + for (i = 0; i < nifnames; i++) {
2.1145 + /* Open input file */
2.1146 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
2.1147 + /*
2.1148 + * Inquire about number of dimensions, variables, global
2.1149 + * attributes and the ID of the unlimited dimension
2.1150 + */
2.1151 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
2.1152 + &input_ngatts, &input_unlimdimid));
2.1153 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
2.1154 + if (strcmp(name, time_name)) {
2.1155 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
2.1156 + exit(4);
2.1157 + }
2.1158 +#ifdef DEBUG
2.1159 + printf("%ld samples in %s\n", (long)len, input_fnames[i]);
2.1160 +#endif /* DEBUG */
2.1161 + wrap_nc(nc_close(input_ncid));
2.1162 + cnt += len;
2.1163 + }
2.1164 + return cnt;
2.1165 +}
2.1166 +
2.1167 +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
2.1168 +{
2.1169 + /*
2.1170 + * Allocate space for arrays of means and standard deviations by
2.1171 + * hour of day.
2.1172 + */
2.1173 + int i;
2.1174 + size_t **cell_samples;
2.1175 + double **means, **stddevs;
2.1176 +
2.1177 + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
2.1178 + perror("alloc_means_stddevs:*meansp");
2.1179 + exit(2);
2.1180 + }
2.1181 + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
2.1182 + perror("alloc_means_stddevs:*stddevsp");
2.1183 + exit(2);
2.1184 + }
2.1185 + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
2.1186 + perror("alloc_means_stddevs:*cell_samplesp");
2.1187 + exit(2);
2.1188 + }
2.1189 +
2.1190 + means = *meansp;
2.1191 + stddevs = *stddevsp;
2.1192 + cell_samples = *cell_samplesp;
2.1193 +
2.1194 + for (i = 0; i < HOURS_PER_DAY; i++) {
2.1195 + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
2.1196 + perror("alloc_means_stddevs:means[i]");
2.1197 + exit(2);
2.1198 + }
2.1199 + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
2.1200 + perror("alloc_means_stddevs:stddevs[i]");
2.1201 + exit(2);
2.1202 + }
2.1203 + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
2.1204 + perror("alloc_means_stddevs:cell_samples[i]");
2.1205 + exit(2);
2.1206 + }
2.1207 + }
2.1208 +
2.1209 + return;
2.1210 +}
2.1211 +
2.1212 +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
2.1213 +{
2.1214 + int hr, i, j, pos;
2.1215 +
2.1216 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
2.1217 + for (i = 0; i < d1; i++) {
2.1218 +#pragma _CRI concurrent
2.1219 + for (j = 0; j < d2; j++) {
2.1220 + pos = i*d2+j;
2.1221 + means[hr][pos] = FillValue;
2.1222 + stddevs[hr][pos] = FillValue;
2.1223 + cell_samples[hr][pos] = 0;
2.1224 + }
2.1225 + }
2.1226 + }
2.1227 +
2.1228 + return;
2.1229 +}
2.1230 +
2.1231 +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
2.1232 +{
2.1233 + int i, j, hr, pos;
2.1234 +
2.1235 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
2.1236 + for (i = 0; i < d1; i++) {
2.1237 +#pragma _CRI concurrent
2.1238 + for (j = 0; j < d2; j++) {
2.1239 + pos = i*d2+j;
2.1240 + cell_samples[hr][pos] = 0;
2.1241 + }
2.1242 + }
2.1243 + }
2.1244 +
2.1245 + return;
2.1246 +}
2.1247 +
2.1248 +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
2.1249 +{
2.1250 + /*
2.1251 + * Free space from arrays of means and standard deviations by
2.1252 + * hour of day.
2.1253 + */
2.1254 + int i;
2.1255 +
2.1256 + for (i = 0; i < HOURS_PER_DAY; i++) {
2.1257 + free(means[i]);
2.1258 + free(stddevs[i]);
2.1259 + free(cell_samples[i]);
2.1260 + }
2.1261 +
2.1262 + free(means);
2.1263 + free(stddevs);
2.1264 + free(cell_samples);
2.1265 +
2.1266 + return;
2.1267 +}
2.1268 +
2.1269 +void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
2.1270 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
2.1271 +{
2.1272 + int i, j, hr, pos;
2.1273 +
2.1274 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
2.1275 +
2.1276 + for (i = 0; i < d1; i++) {
2.1277 +#pragma _CRI concurrent
2.1278 + for (j = 0; j < d2; j++) {
2.1279 + pos = i*d2+j;
2.1280 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
2.1281 + if (cell_samples[hr][pos] == 0)
2.1282 + means[hr][pos] = (double)val[pos];
2.1283 + else
2.1284 + means[hr][pos] += (double)val[pos];
2.1285 + ++cell_samples[hr][pos];
2.1286 + }
2.1287 + }
2.1288 + }
2.1289 +
2.1290 + return;
2.1291 +}
2.1292 +
2.1293 +void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
2.1294 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
2.1295 +{
2.1296 + int i, j, hr, pos;
2.1297 +
2.1298 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
2.1299 +
2.1300 + for (i = 0; i < d1; i++) {
2.1301 +#pragma _CRI concurrent
2.1302 + for (j = 0; j < d2; j++) {
2.1303 + pos = i*d2+j;
2.1304 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
2.1305 + if (cell_samples[hr][pos] == 0)
2.1306 + means[hr][pos] = val[pos];
2.1307 + else
2.1308 + means[hr][pos] += val[pos];
2.1309 + ++cell_samples[hr][pos];
2.1310 + }
2.1311 + }
2.1312 + }
2.1313 +
2.1314 + return;
2.1315 +}
2.1316 +
2.1317 +
2.1318 +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
2.1319 +{
2.1320 + int i, j, hr, pos;
2.1321 +
2.1322 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
2.1323 + for (i = 0; i < d1; i++) {
2.1324 +#pragma _CRI concurrent
2.1325 + for (j = 0; j < d2; j++) {
2.1326 + pos = i*d2+j;
2.1327 + if (cell_samples[hr][pos] != 0) {
2.1328 + means[hr][pos] /= (double)cell_samples[hr][pos];
2.1329 + }
2.1330 + }
2.1331 + }
2.1332 + }
2.1333 +
2.1334 + return;
2.1335 +}
2.1336 +
2.1337 +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
2.1338 + char FillFlag, float FillValue, double **means, double **stddevs,
2.1339 + size_t **cell_samples)
2.1340 +{
2.1341 + int i, j, hr, pos;
2.1342 + double delta;
2.1343 +
2.1344 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
2.1345 +
2.1346 + for (i = 0; i < d1; i++) {
2.1347 +#pragma _CRI concurrent
2.1348 + for (j = 0; j < d2; j++) {
2.1349 + pos = i*d2+j;
2.1350 + if (!FillFlag || (FillFlag && val[pos] != FillValue
2.1351 + && means[hr][pos] != FillValue)) {
2.1352 + delta = means[hr][pos] - (double)val[pos];
2.1353 + if (cell_samples[hr][pos] == 0)
2.1354 + stddevs[hr][pos] = delta * delta;
2.1355 + else
2.1356 + stddevs[hr][pos] += delta * delta;
2.1357 + ++cell_samples[hr][pos];
2.1358 + }
2.1359 + }
2.1360 + }
2.1361 +
2.1362 + return;
2.1363 +}
2.1364 +
2.1365 +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
2.1366 + char FillFlag, float FillValue, double **means, double **stddevs,
2.1367 + size_t **cell_samples)
2.1368 +{
2.1369 + int i, j, hr, pos;
2.1370 + double delta;
2.1371 +
2.1372 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
2.1373 +
2.1374 + for (i = 0; i < d1; i++) {
2.1375 +#pragma _CRI concurrent
2.1376 + for (j = 0; j < d2; j++) {
2.1377 + pos = i*d2+j;
2.1378 + if (!FillFlag || (FillFlag && val[pos] != FillValue
2.1379 + && means[hr][pos] != FillValue)) {
2.1380 + delta = means[hr][pos] - val[pos];
2.1381 + if (cell_samples[hr][pos] == 0)
2.1382 + stddevs[hr][pos] = delta * delta;
2.1383 + else
2.1384 + stddevs[hr][pos] += delta * delta;
2.1385 + ++cell_samples[hr][pos];
2.1386 + }
2.1387 + }
2.1388 + }
2.1389 +
2.1390 + return;
2.1391 +}
2.1392 +
2.1393 +
2.1394 +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
2.1395 +{
2.1396 + int i, j, hr, pos;
2.1397 +
2.1398 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
2.1399 + for (i = 0; i < d1; i++) {
2.1400 +#pragma _CRI concurrent
2.1401 + for (j = 0; j < d2; j++) {
2.1402 + pos = i*d2+j;
2.1403 + if (cell_samples[hr][pos] != 0) {
2.1404 + stddevs[hr][pos] /= (double)cell_samples[hr][pos];
2.1405 + stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
2.1406 + }
2.1407 + }
2.1408 + }
2.1409 + }
2.1410 +
2.1411 + return;
2.1412 +}
2.1413 +
2.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)
2.1415 +{
2.1416 + int i, ts;
2.1417 + int varid, *mcsec;
2.1418 + size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.1419 + char name[NC_MAX_NAME+1];
2.1420 + void *val;
2.1421 +
2.1422 + /* Allocate space for timeslice */
2.1423 + val = alloc_var(nctype, (d1 * d2));
2.1424 +
2.1425 + for (i = 0; i < nifnames; i++) {
2.1426 +#ifdef DEBUG
2.1427 + printf("\tOpening %s", input_fnames[i]);
2.1428 + if (ndims > 3) printf(" to retrieve level %d", level);
2.1429 + if (stat_type == MEAN_TYPE)
2.1430 + printf(" for computing means\n");
2.1431 + else
2.1432 + printf(" for computing stddevs\n");
2.1433 +#endif /* DEBUG */
2.1434 + /* Open input file */
2.1435 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
2.1436 + /*
2.1437 + * Inquire about number of dimensions, variables, global
2.1438 + * attributes and the ID of the unlimited dimension
2.1439 + */
2.1440 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
2.1441 + &input_ngatts, &input_unlimdimid));
2.1442 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
2.1443 + if (strcmp(name, time_name)) {
2.1444 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
2.1445 + exit(4);
2.1446 + }
2.1447 + /* Get seconds of day variable */
2.1448 + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
2.1449 + if (!(mcsec = (int *)malloc(len * sizeof(int)))) {
2.1450 + perror("read_and_compute:mcsec");
2.1451 + exit(2);
2.1452 + }
2.1453 + wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));
2.1454 + /* Get variable ID for requested variable */
2.1455 + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
2.1456 + /* Retrieve time slice of desired variable */
2.1457 + for (ts = 0; ts < len; ts++) {
2.1458 + start[0] = ts;
2.1459 + count[0] = 1;
2.1460 + start[(ndims-2)] = 0;
2.1461 + count[(ndims-2)] = d1;
2.1462 + start[(ndims-1)] = 0;
2.1463 + count[(ndims-1)] = d2;
2.1464 + if (ndims == 4) {
2.1465 + start[1] = level;
2.1466 + count[1] = 1;
2.1467 + }
2.1468 + switch(nctype) {
2.1469 + case NC_FLOAT:
2.1470 + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));
2.1471 + if (stat_type == MEAN_TYPE)
2.1472 + add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
2.1473 + else
2.1474 + add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
2.1475 + break;
2.1476 + case NC_DOUBLE:
2.1477 + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));
2.1478 + if (stat_type == MEAN_TYPE)
2.1479 + add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
2.1480 + else
2.1481 + add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
2.1482 + break;
2.1483 + default:
2.1484 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.1485 + exit(3);
2.1486 + }
2.1487 + }
2.1488 +
2.1489 + /* Free mcsec vector */
2.1490 + free(mcsec);
2.1491 + /* Close input file */
2.1492 + wrap_nc(nc_close(input_ncid));
2.1493 + }
2.1494 + /* Divide sums by number of samples */
2.1495 + if (stat_type == MEAN_TYPE)
2.1496 + divide_means(d1, d2, means, cell_samples);
2.1497 + else
2.1498 + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
2.1499 +
2.1500 + /* Free working variable array */
2.1501 + free(val);
2.1502 +
2.1503 + return;
2.1504 +}
2.1505 +
2.1506 +float *double_to_float(size_t len, double *dvar)
2.1507 +{
2.1508 + int i;
2.1509 + float *fvar;
2.1510 +
2.1511 + if (!(fvar = (float *)malloc(len * sizeof(float)))) {
2.1512 + perror("double_to_float:fvar");
2.1513 + exit(2);
2.1514 + }
2.1515 +
2.1516 + for (i = 0; i < len; i++)
2.1517 + fvar[i] = (float)dvar[i];
2.1518 +
2.1519 + return fvar;
2.1520 +}
2.1521 +
2.1522 +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
2.1523 + int level, int ndims, double **var_hours)
2.1524 +{
2.1525 + /* Output dimensions should be one larger than input dimensions */
2.1526 + int i, hr;
2.1527 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.1528 + float *fvar = NULL;
2.1529 +
2.1530 + if (nctype == NC_FLOAT) {
2.1531 + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
2.1532 + perror("write_var_hours:fvar");
2.1533 + exit(2);
2.1534 + }
2.1535 + }
2.1536 +
2.1537 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
2.1538 + start[0] = hr;
2.1539 + count[0] = 1; /* time */
2.1540 + start[(ndims-2)] = 0;
2.1541 + count[(ndims-2)] = d1;
2.1542 + start[(ndims-1)] = 0;
2.1543 + count[(ndims-1)] = d2;
2.1544 + if (ndims == 4) {
2.1545 + start[1] = level;
2.1546 + count[1] = 1;
2.1547 + }
2.1548 + switch (nctype) {
2.1549 + case NC_FLOAT:
2.1550 + for (i = 0; i < (d1 * d2); i++)
2.1551 + fvar[i] = (float)var_hours[hr][i];
2.1552 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
2.1553 + break;
2.1554 + case NC_DOUBLE:
2.1555 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
2.1556 + break;
2.1557 + default:
2.1558 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
2.1559 + exit(3);
2.1560 + }
2.1561 + }
2.1562 +
2.1563 + if (nctype == NC_FLOAT)
2.1564 + free(fvar);
2.1565 +
2.1566 + return;
2.1567 +}
2.1568 +
2.1569 +void compute_stats(int nifnames, char **input_fnames)
2.1570 +{
2.1571 + int j, k, nlevels;
2.1572 + size_t d1, d2, **cell_samples;
2.1573 + double **means, **stddevs;
2.1574 + struct var *in_vnode, *out_vnode;
2.1575 +
2.1576 + if (input_var_head) {
2.1577 + in_vnode = input_var_head;
2.1578 + /* March through non-metadata variables to compute stats */
2.1579 + for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
2.1580 + if (!in_vnode->metadata) {
2.1581 + /* Make sure time is really the first dimension */
2.1582 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
2.1583 + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
2.1584 + exit(3);
2.1585 + }
2.1586 + /* Figure out input dimensions */
2.1587 + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
2.1588 + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
2.1589 + exit(3);
2.1590 + }
2.1591 + /* Assume 2D output; find dimensions */
2.1592 + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
2.1593 + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
2.1594 + if (in_vnode->ndims == 3)
2.1595 + nlevels = 1;
2.1596 + else
2.1597 + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
2.1598 + /* Allocate working space for means and stddevs */
2.1599 + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
2.1600 + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
2.1601 + printf("Computing statistics for %s\n",
2.1602 + in_vnode->name);
2.1603 + /* Compute means and stddevs, then write them */
2.1604 + for (k = 0; k < nlevels; k++) {
2.1605 + /* Read and compute means */
2.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);
2.1607 + /* Find corresponding output variable on the mean output file */
2.1608 + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
2.1609 + /* Write out the means for this variable */
2.1610 + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
2.1611 + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
2.1612 + reset_cell_samples(d1, d2, cell_samples);
2.1613 + /* Read and compute stddevs using means */
2.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);
2.1615 + /* Find corresponding output variable on the stddev output file */
2.1616 + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
2.1617 + /* Write out stddevs for this variable */
2.1618 + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
2.1619 + }
2.1620 +
2.1621 + /* Free working space for means and stddevs */
2.1622 + free_means_stddevs(means, stddevs, cell_samples);
2.1623 + }
2.1624 + in_vnode = in_vnode->next;
2.1625 + }
2.1626 + }
2.1627 + return;
2.1628 +}
2.1629 +
2.1630 +void usage(char *program)
2.1631 +{
2.1632 + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
2.1633 + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
2.1634 + exit(4);
2.1635 +}
2.1636 +
2.1637 +int main(int argc, char **argv)
2.1638 +{
2.1639 + int i, ic, nifnames;
2.1640 + size_t len, pos, nsamples;
2.1641 + char *mfname, *sfname, **ifnames, *flist;
2.1642 +
2.1643 + ifnames = NULL;
2.1644 + mfname = sfname = flist = NULL;
2.1645 + nifnames = 0;
2.1646 +
2.1647 +#ifdef DEBUG
2.1648 + printf("Number of metadata variables (nmvars) = %d\n", nmvars);
2.1649 + fflush(stdout);
2.1650 +#endif /* DEBUG */
2.1651 +
2.1652 + /* check command-line flags and switches */
2.1653 + while ((ic = getopt(argc, argv, "m:s:")) != -1) {
2.1654 + switch(ic) {
2.1655 + case 'm':
2.1656 + if (!(mfname = strdup(optarg))) {
2.1657 + perror("mfname");
2.1658 + exit(2);
2.1659 + }
2.1660 + break;
2.1661 + case 's':
2.1662 + if (!(sfname = strdup(optarg))) {
2.1663 + perror("sfname");
2.1664 + exit(2);
2.1665 + }
2.1666 + break;
2.1667 + }
2.1668 + }
2.1669 +
2.1670 + if (!mfname) {
2.1671 + fprintf(stderr, "Output file name for writing means required.\n");
2.1672 + usage(argv[0]);
2.1673 + }
2.1674 + if (!sfname) {
2.1675 + fprintf(stderr, "Output file name for writing standard deviations required.\n");
2.1676 + usage(argv[0]);
2.1677 + }
2.1678 +
2.1679 + if (optind < argc) {
2.1680 + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
2.1681 + perror("ifnames");
2.1682 + exit(2);
2.1683 + }
2.1684 + for (i = optind; i < argc; i++) {
2.1685 + if (!(ifnames[nifnames++] = strdup(argv[i]))) {
2.1686 + perror("ifnames[nifnames++]");
2.1687 + exit(2);
2.1688 + }
2.1689 + }
2.1690 + }
2.1691 + else {
2.1692 + fprintf(stderr, "At least one input file name is required.\n");
2.1693 + usage(argv[0]);
2.1694 + }
2.1695 +
2.1696 +
2.1697 + /*
2.1698 + * Build list of input files to be included in the global history
2.1699 + * attribute on the two outputs files.
2.1700 + */
2.1701 + for (i = len = 0; i < nifnames; i++)
2.1702 + len += strlen(ifnames[i]);
2.1703 + len += nifnames + 1; /* space in front of every name + null terminator */
2.1704 + if (!(flist = (char *)malloc(len * sizeof(char)))) {
2.1705 + perror("flist");
2.1706 + exit(2);
2.1707 + }
2.1708 + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
2.1709 + sprintf(&flist[pos], " %s", ifnames[i]);
2.1710 +#ifdef DEBUG
2.1711 + printf("flist=%s\n", flist);
2.1712 + fflush(stdout);
2.1713 +#endif /* DEBUG */
2.1714 +
2.1715 + /* Open every file to count up the number of time samples. */
2.1716 + nsamples = count_samples(nifnames, ifnames);
2.1717 +#ifdef DEBUG
2.1718 + printf("Number of samples across all files = %ld\n", (long)nsamples);
2.1719 + fflush(stdout);
2.1720 +#endif /* DEBUG */
2.1721 +
2.1722 + /*
2.1723 + * Open the *last* input file and the two output files (for mean and
2.1724 + * standard deviation). Define dimensions, variables, and attributes
2.1725 + * for the two output files based on the *last* input files. The
2.1726 + * last file is used because some metadata variables will contain
2.1727 + * only values from the last time slice from the period over which
2.1728 + * the statistics are computed.
2.1729 + */
2.1730 + open_inout(ifnames[0], ifnames[(nifnames-1)], mfname, sfname, flist,
2.1731 + nsamples);
2.1732 +
2.1733 + compute_stats(nifnames, ifnames);
2.1734 +
2.1735 + /* Close the two output files */
2.1736 + wrap_nc(nc_close(mean_ncid));
2.1737 + wrap_nc(nc_close(stddev_ncid));
2.1738 +
2.1739 + return 0;
2.1740 +}