1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/Makefile Wed Sep 26 17:16:40 2007 -0400
1.3 @@ -0,0 +1,38 @@
1.4 +CC=gcc
1.5 +#
1.6 +# robin1 and phoenix
1.7 +LIBS=-L$(NETCDF)/lib -lnetcdf -lm
1.8 +# Penguins
1.9 +#LIBS=-L/usr/local/netcdf/netcdf-3.6.1-gcc+pgi/lib -lnetcdf -lm
1.10 +#
1.11 +# robin1 and phoenix
1.12 +CPPFLAGS=-I$(NETCDF)/include
1.13 +# Penguins
1.14 +#CPPFLAGS=-I/usr/local/netcdf/netcdf-3.6.1-gcc+pgi/include
1.15 +# phoenix
1.16 +#CC=cc
1.17 +#CFLAGS=-O -h list=m $(CPPFLAGS)
1.18 +# robin1 and Penguins
1.19 +CFLAGS=-g -Wall -O $(CPPFLAGS)
1.20 +
1.21 +all: h1_summary h1_summary2 add_total_fields
1.22 +
1.23 +h1_summary: h1_summary.o
1.24 + $(CC) $(CFLAGS) -o $@ h1_summary.o $(LIBS)
1.25 +
1.26 +h1_summary2: h1_summary2.o
1.27 + $(CC) $(CFLAGS) -o $@ h1_summary2.o $(LIBS)
1.28 +
1.29 +add_total_fields: add_total_fields.o
1.30 + $(CC) $(CFLAGS) -o $@ add_total_fields.o $(LIBS)
1.31 +
1.32 +clean:
1.33 + $(RM) -f h1_summary.o h1_summary
1.34 + $(RM) -f h1_summary2.o h1_summary2
1.35 + $(RM) -f add_total_fields.o add_total_fields
1.36 +
1.37 +install: all
1.38 + cp -p h1_summary $(HOME)/bin/
1.39 + cp -p h1_summary2 $(HOME)/bin/
1.40 + cp -p add_total_fields $(HOME)/bin/
1.41 +
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/add_total_fields.c Wed Sep 26 17:16:40 2007 -0400
2.3 @@ -0,0 +1,548 @@
2.4 +#include <stdio.h>
2.5 +#include <stdlib.h>
2.6 +#include <string.h>
2.7 +#include <math.h>
2.8 +#include "netcdf.h"
2.9 +
2.10 +/*
2.11 + * Loop through all history tapes from the Community Land Model adding new
2.12 + * fields that are totals of the multi-level fields included in total_vars[]
2.13 + * below.
2.14 + */
2.15 +
2.16 +static char *total_vars[] = {
2.17 + "SOILICE",
2.18 + "SOILLIQ"
2.19 +};
2.20 +static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars);
2.21 +static char *combo_vars[] = {
2.22 + "LATENT", /* FCTR + FCEV + FGEV */
2.23 + "ET", /* QVEGE + QVEGT + QSOIL */
2.24 + "ALL_SKY_ALBEDO", /* FSR / FSDS */
2.25 + "BLACK_SKY_ALBEDO", /* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */
2.26 + "NETRAD", /* FSA - FIRA */
2.27 +};
2.28 +static char *combo_long_names[] = {
2.29 + "latent heat flux",
2.30 + "evapotranspiration",
2.31 + "surface all-sky albedo",
2.32 + "surface black-sky albedo",
2.33 + "net radiation"
2.34 +};
2.35 +static char *combo_units[] = {
2.36 + "watt/m^2",
2.37 + "m/s",
2.38 + "unitless",
2.39 + "unitless",
2.40 + "watt/m^2"
2.41 +};
2.42 +static float combo_FillValues[] = {
2.43 + 1.e+36,
2.44 + 1.e+36,
2.45 + 1.e+36,
2.46 + 1.e+36,
2.47 + 1.e+36
2.48 +};
2.49 +static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars);
2.50 +static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total ";
2.51 +static char *time_name = "time", *lon_name = "lon", *lat_name = "lat";
2.52 +static char *long_name_name = "long_name", *units_name = "units",
2.53 + *cell_method_name = "cell_method", *FillValue_name = "_FillValue",
2.54 + *missing_value_name = "missing_value";
2.55 +static char *cell_method = "time: mean";
2.56 +
2.57 +void wrap_nc(int status)
2.58 +{
2.59 + if (status != NC_NOERR) {
2.60 + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
2.61 + exit(1);
2.62 + }
2.63 +
2.64 + return;
2.65 +}
2.66 +
2.67 +void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1,
2.68 + size_t d2, char FillFlag, float FillValue, size_t tslice, float *var,
2.69 + double *tvar)
2.70 +{
2.71 + int i, j;
2.72 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.73 +
2.74 + /* Initialize tvar to FillValue */
2.75 + for (i = 0; i < (d1 * d2); i++)
2.76 + tvar[i] = (double)FillValue;
2.77 +
2.78 +
2.79 + start[0] = tslice;
2.80 + count[0] = 1;
2.81 + start[2] = 0;
2.82 + count[2] = d1;
2.83 + start[3] = 0;
2.84 + count[3] = d2;
2.85 +
2.86 + /* March through the levels, totalling as we go */
2.87 + for (j = 0; j < nlevels; j++) {
2.88 + start[1] = j;
2.89 + count[1] = 1;
2.90 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
2.91 + for (i = 0; i < (d1 * d2); i++) {
2.92 + if (tvar[i] == (double)FillValue)
2.93 + tvar[i] = (double)var[i];
2.94 + else
2.95 + tvar[i] += (double)var[i];
2.96 + }
2.97 + }
2.98 +
2.99 + return;
2.100 +}
2.101 +
2.102 +void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)
2.103 +{
2.104 + int i;
2.105 + float **svars;
2.106 +
2.107 + if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {
2.108 + perror("alloc_svars:svarsp");
2.109 + exit(2);
2.110 + }
2.111 + svars = *svarsp;
2.112 + for (i = 0; i < nsvars; i++) {
2.113 + if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) {
2.114 + perror("alloc_svars:svars[i]");
2.115 + exit(2);
2.116 + }
2.117 + }
2.118 + return;
2.119 +}
2.120 +
2.121 +void free_svars(size_t nsvars, float **svars)
2.122 +{
2.123 + int i;
2.124 +
2.125 + for (i = 0; i < nsvars; i++)
2.126 + free(svars[i]);
2.127 + free(svars);
2.128 + return;
2.129 +}
2.130 +
2.131 +void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,
2.132 + float *var)
2.133 +{
2.134 + int varid;
2.135 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.136 +
2.137 + start[0] = tslice;
2.138 + count[0] = 1;
2.139 + start[1] = 0;
2.140 + count[1] = d1;
2.141 + start[2] = 0;
2.142 + count[2] = d2;
2.143 +
2.144 + /* Get variable id */
2.145 + wrap_nc(nc_inq_varid(ncid, name, &varid));
2.146 + /* Assume the correct size and read the variable */
2.147 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
2.148 +
2.149 + return;
2.150 +}
2.151 +
2.152 +void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,
2.153 + int num, double *tvar)
2.154 +{
2.155 + int i;
2.156 + size_t nsvars;
2.157 + float **svars; /* source variables */
2.158 + double denom;
2.159 +
2.160 + /* Initialize tvar to FillValue */
2.161 + for (i = 0; i < (d1 * d2); i++)
2.162 + tvar[i] = (double)combo_FillValues[num];
2.163 +
2.164 + switch (num) {
2.165 + case 0: /* LATENT */
2.166 + nsvars = 3;
2.167 + alloc_svars(nsvars, d1, d2, &svars);
2.168 + /* read timeslices */
2.169 + read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]);
2.170 + read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]);
2.171 + read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]);
2.172 + /* compute new variable */
2.173 + for (i = 0; i < (d1 * d2); i++) {
2.174 + if (svars[0][i] != combo_FillValues[num] &&
2.175 + svars[1][i] != combo_FillValues[num] &&
2.176 + svars[2][i] != combo_FillValues[num]) {
2.177 + tvar[i] = (double)svars[0][i]
2.178 + + (double)svars[1][i]
2.179 + + (double)svars[2][i];
2.180 + }
2.181 + }
2.182 + free_svars(nsvars, svars);
2.183 + break;
2.184 + case 1: /* ET */
2.185 + nsvars = 3;
2.186 + alloc_svars(nsvars, d1, d2, &svars);
2.187 + /* read timeslices */
2.188 + read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]);
2.189 + read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]);
2.190 + read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]);
2.191 + /* compute new variable */
2.192 + for (i = 0; i < (d1 * d2); i++) {
2.193 + if (svars[0][i] != combo_FillValues[num] &&
2.194 + svars[1][i] != combo_FillValues[num] &&
2.195 + svars[2][i] != combo_FillValues[num]) {
2.196 + tvar[i] = (double)svars[0][i]
2.197 + + (double)svars[1][i]
2.198 + + (double)svars[2][i];
2.199 + }
2.200 + }
2.201 + free_svars(nsvars, svars);
2.202 + break;
2.203 + case 2: /* ALL_SKY_ALBEDO */
2.204 + nsvars = 2;
2.205 + alloc_svars(nsvars, d1, d2, &svars);
2.206 + /* read timeslices */
2.207 + read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]);
2.208 + read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]);
2.209 + /* compute new variable */
2.210 + for (i = 0; i < (d1 * d2); i++) {
2.211 + if (svars[0][i] != combo_FillValues[num] &&
2.212 + svars[1][i] != combo_FillValues[num] &&
2.213 + (svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) {
2.214 + tvar[i] = (double)svars[0][i]
2.215 + / (double)svars[1][i];
2.216 + }
2.217 + }
2.218 + free_svars(nsvars, svars);
2.219 + break;
2.220 + case 3: /* BLACK_SKY_ALBEDO */
2.221 + nsvars = 4;
2.222 + alloc_svars(nsvars, d1, d2, &svars);
2.223 + /* read timeslices */
2.224 + read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]);
2.225 + read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]);
2.226 + read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]);
2.227 + read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]);
2.228 + /* compute new variable */
2.229 + for (i = 0; i < (d1 * d2); i++) {
2.230 + if (svars[0][i] != combo_FillValues[num] &&
2.231 + svars[1][i] != combo_FillValues[num] &&
2.232 + svars[2][i] != combo_FillValues[num] &&
2.233 + svars[3][i] != combo_FillValues[num]) {
2.234 + denom = (double)svars[2][i] + (double)svars[3][i];
2.235 + if (denom > 1.e-35 || denom < -1.e-35) {
2.236 + tvar[i] = ((double)svars[0][i]
2.237 + + (double)svars[1][i]) /
2.238 + ((double)svars[2][i]
2.239 + + (double)svars[3][i]);
2.240 + }
2.241 + }
2.242 + }
2.243 + free_svars(nsvars, svars);
2.244 + break;
2.245 + case 4: /* NETRAD */
2.246 + nsvars = 2;
2.247 + alloc_svars(nsvars, d1, d2, &svars);
2.248 + /* read timeslices */
2.249 + read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]);
2.250 + read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]);
2.251 + /* compute new variable */
2.252 + for (i = 0; i < (d1 * d2); i++) {
2.253 + if (svars[0][i] != combo_FillValues[num] &&
2.254 + svars[1][i] != combo_FillValues[num]) {
2.255 + tvar[i] = (double)svars[0][i]
2.256 + - (double)svars[1][i];
2.257 + }
2.258 + }
2.259 + free_svars(nsvars, svars);
2.260 + break;
2.261 + default:
2.262 + fprintf(stderr, "Unknown combination variable %d!\n", num);
2.263 + exit(3);
2.264 + }
2.265 +
2.266 + return;
2.267 +}
2.268 +
2.269 +void write_timeslice(int ncid, int varid, size_t d1, size_t d2,
2.270 + size_t tslice, float *var)
2.271 +{
2.272 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
2.273 +
2.274 + start[0] = tslice;
2.275 + count[0] = 1;
2.276 + start[1] = 0;
2.277 + count[1] = d1;
2.278 + start[2] = 0;
2.279 + count[2] = d2;
2.280 +
2.281 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));
2.282 +
2.283 + return;
2.284 +}
2.285 +
2.286 +void usage(char *program)
2.287 +{
2.288 + fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",
2.289 + program);
2.290 + fprintf(stderr, "Requires at least one history output file name\n");
2.291 + exit(5);
2.292 +}
2.293 +
2.294 +int main(int argc, char **argv)
2.295 +{
2.296 + char **ifnames;
2.297 + int i, j, k, nifnames;
2.298 + int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts,
2.299 + dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS],
2.300 + new_varid, lat_dimid, lon_dimid;
2.301 + nc_type xtype, atype;
2.302 + char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1],
2.303 + att_name[NC_MAX_NAME+1], *longname, *new_longname;
2.304 + size_t tlen, ts, alen, nlevels, d1, d2;
2.305 + float FillValue;
2.306 + char FillFlag;
2.307 + float *var;
2.308 + double *tvar;
2.309 +
2.310 +#ifdef DEBUG
2.311 + printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);
2.312 + printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);
2.313 +#endif /* DEBUG */
2.314 +
2.315 + /* Require at least one parameter (the file on which to work) */
2.316 + if (argc < 2)
2.317 + usage(argv[0]);
2.318 +
2.319 + /* Put the list of filenames in an array of strings */
2.320 + if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {
2.321 + perror("ifnames");
2.322 + exit(2);
2.323 + }
2.324 + nifnames = 0;
2.325 + for (i = 1; i < argc; i++) {
2.326 + if (!(ifnames[nifnames++] = strdup(argv[i]))) {
2.327 + perror("ifnames[nifnames++]");
2.328 + exit(2);
2.329 + }
2.330 + }
2.331 +
2.332 + /*
2.333 + * March through all input files adding new variables for totals
2.334 + * across levels.
2.335 + */
2.336 + for (i = 0; i < nifnames; i++) {
2.337 + printf("Processing %s\n", ifnames[i]);
2.338 + /* Open file */
2.339 + wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));
2.340 + /*
2.341 + * Inquire about number of dimensions, variables, global
2.342 + * attributes and the ID of the unlimited dimension
2.343 + */
2.344 + wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid));
2.345 + wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen));
2.346 + if (strcmp(name, time_name)) {
2.347 + fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]);
2.348 + exit(4);
2.349 + }
2.350 + /* Retrieve lat and lon dimension IDs */
2.351 + wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid));
2.352 + wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid));
2.353 + /* Put file into define mode */
2.354 + wrap_nc(nc_redef(ncid));
2.355 + /*
2.356 + * Define all of the new variables first to improve
2.357 + * performance. This means some calls are made twice.
2.358 + */
2.359 + /* First, add the new total variables */
2.360 + for (j = 0; j < nmtotal_vars; j++) {
2.361 + /* Figure out source variable information */
2.362 + wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
2.363 + wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
2.364 + /* Make sure it has the correct number of dimensions */
2.365 + if (ndims != 4) {
2.366 + fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);
2.367 + exit(4);
2.368 + }
2.369 + /* Make sure time is the first dimension */
2.370 + if (dimids[0] != unlimdimid) {
2.371 + fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name);
2.372 + exit(4);
2.373 + }
2.374 + /* Make sure it is a float type */
2.375 + if (xtype != NC_FLOAT) {
2.376 + fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");
2.377 + exit(4);
2.378 + }
2.379 + /* Set dimensions for new variable */
2.380 + new_ndims = 3;
2.381 + new_dimids[0] = unlimdimid;
2.382 + new_dimids[1] = dimids[2];
2.383 + new_dimids[2] = dimids[3];
2.384 + /* Define new variable */
2.385 + sprintf(new_name, "%s%s", varname_total_prefix, name);
2.386 + printf("\tAdding variable %s\n", new_name);
2.387 + wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims,
2.388 + new_dimids, &new_varid));
2.389 + /* Copy attributes from source variable to new one */
2.390 + for (k = 0; k < natts; k++) {
2.391 + wrap_nc(nc_inq_attname(ncid, varid, k, att_name));
2.392 + if (!strcmp(att_name, "long_name")) {
2.393 + wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen));
2.394 + if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) {
2.395 + perror("longname");
2.396 + exit(2);
2.397 + }
2.398 + wrap_nc(nc_get_att_text(ncid, varid, att_name, longname));
2.399 + longname[alen] = '\0';
2.400 + if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) {
2.401 + perror("new_longname");
2.402 + exit(2);
2.403 + }
2.404 + sprintf(new_longname, "%s%s", longname_total_prefix, longname);
2.405 + wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname));
2.406 + free(new_longname);
2.407 + free(longname);
2.408 + }
2.409 + else
2.410 + wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));
2.411 + }
2.412 + }
2.413 + /* Second, add the new total variables */
2.414 + for (j = 0; j < nmcombo_vars; j++) {
2.415 + /* Set dimensions for new variable */
2.416 + new_ndims = 3;
2.417 + new_dimids[0] = unlimdimid;
2.418 + new_dimids[1] = lat_dimid;
2.419 + new_dimids[2] = lon_dimid;
2.420 + /* Define new variable */
2.421 + printf("\tAdding variable %s\n", combo_vars[j]);
2.422 + wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT,
2.423 + new_ndims, new_dimids, &new_varid));
2.424 + /* Add attributes */
2.425 + wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name,
2.426 + strlen(combo_long_names[j]),
2.427 + combo_long_names[j]));
2.428 + wrap_nc(nc_put_att_text(ncid, new_varid, units_name,
2.429 + strlen(combo_units[j]), combo_units[j]));
2.430 + wrap_nc(nc_put_att_text(ncid, new_varid,
2.431 + cell_method_name, strlen(cell_method),
2.432 + cell_method));
2.433 + wrap_nc(nc_put_att_float(ncid, new_varid,
2.434 + FillValue_name, NC_FLOAT, 1,
2.435 + &combo_FillValues[j]));
2.436 + wrap_nc(nc_put_att_float(ncid, new_varid,
2.437 + missing_value_name, NC_FLOAT, 1,
2.438 + &combo_FillValues[j]));
2.439 + }
2.440 + /* Leave define mode */
2.441 + wrap_nc(nc_enddef(ncid));
2.442 + /* Now actually compute and write out the new total variables */
2.443 + for (j = 0; j < nmtotal_vars; j++) {
2.444 + /* Figure out source variable information */
2.445 + wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
2.446 + wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
2.447 + /* Check for _FillValue */
2.448 + FillFlag = 0;
2.449 + FillValue = 0.;
2.450 + if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) {
2.451 + if (atype == NC_FLOAT && alen) {
2.452 + wrap_nc(nc_get_att_float(ncid, varid,
2.453 + FillValue_name, &FillValue));
2.454 + FillFlag = 1;
2.455 + }
2.456 + }
2.457 + /* Set dimensions for new variable */
2.458 + new_ndims = 3;
2.459 + new_dimids[0] = unlimdimid;
2.460 + new_dimids[1] = dimids[2];
2.461 + new_dimids[2] = dimids[3];
2.462 +
2.463 + sprintf(new_name, "%s%s", varname_total_prefix, name);
2.464 + printf("\tComputing and writing %s\n", new_name);
2.465 + /* Retrieve the new varid again */
2.466 + wrap_nc(nc_inq_varid(ncid, new_name, &new_varid));
2.467 + /*
2.468 + * Figure out dimensions and total size
2.469 + * for a single time slice
2.470 + */
2.471 + wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels));
2.472 + wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1));
2.473 + wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2));
2.474 +
2.475 + /*
2.476 + * Allocate space for reading in time slice and for
2.477 + * the sum
2.478 + */
2.479 + if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
2.480 + perror("var");
2.481 + exit(2);
2.482 + }
2.483 + if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
2.484 + perror("tvar");
2.485 + exit(2);
2.486 + }
2.487 +
2.488 + /* Read timeslices, total, and write out new timeslices */
2.489 + for (ts = 0; ts < tlen; ts++) {
2.490 + read_total_timeslice(ncid, varid,
2.491 + nlevels, d1, d2, FillFlag, FillValue,
2.492 + ts, var, tvar);
2.493 + for (k = 0; k < (d1*d2); k++)
2.494 + var[k] = (float)tvar[k];
2.495 + write_timeslice(ncid, new_varid, d1, d2,
2.496 + ts, var);
2.497 + }
2.498 +
2.499 + free(var);
2.500 + free(tvar);
2.501 + }
2.502 + /* Now actually compute and write out the new combo variables */
2.503 + for (j = 0; j < nmcombo_vars; j++) {
2.504 + /* Set dimensions for new variable */
2.505 + new_ndims = 3;
2.506 + new_dimids[0] = unlimdimid;
2.507 + new_dimids[1] = lat_dimid;
2.508 + new_dimids[2] = lon_dimid;
2.509 +
2.510 + printf("\tComputing and writing %s\n", combo_vars[j]);
2.511 + /* Retrieve the new varid again */
2.512 + wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid));
2.513 + /*
2.514 + * Figure out dimensions and total size
2.515 + * for a single time slice
2.516 + */
2.517 + wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));
2.518 + wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));
2.519 +
2.520 + /*
2.521 + * Allocate space for reading in time slice and for
2.522 + * the sum
2.523 + */
2.524 + if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
2.525 + perror("var");
2.526 + exit(2);
2.527 + }
2.528 + if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
2.529 + perror("tvar");
2.530 + exit(2);
2.531 + }
2.532 +
2.533 + /* Read timeslices, compute new variable, and write */
2.534 + for (ts = 0; ts < tlen; ts++) {
2.535 + read_compute_timeslice(ncid, d1, d2, ts,
2.536 + j, tvar);
2.537 + for (k = 0; k < (d1*d2); k++)
2.538 + var[k] = (float)tvar[k];
2.539 + write_timeslice(ncid, new_varid, d1, d2,
2.540 + ts, var);
2.541 + }
2.542 +
2.543 + free(var);
2.544 + free(tvar);
2.545 + }
2.546 + /* Close file */
2.547 + wrap_nc(nc_close(ncid));
2.548 + }
2.549 +
2.550 + return 0;
2.551 +}
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/h1_summary.c Wed Sep 26 17:16:40 2007 -0400
3.3 @@ -0,0 +1,1525 @@
3.4 +#include <stdio.h>
3.5 +#include <unistd.h>
3.6 +#include <stdlib.h>
3.7 +#include <string.h>
3.8 +#include <math.h>
3.9 +#include "netcdf.h"
3.10 +
3.11 +/*
3.12 + * Loop through one month of hourly history tapes from the Community Land Model
3.13 + * and generate monthly statistics (means and standard deviations) of fields
3.14 + * by hour of day.
3.15 + */
3.16 +
3.17 +#define MEAN_TYPE 0
3.18 +#define STDDEV_TYPE 1
3.19 +
3.20 +#define HOURS_PER_DAY 24
3.21 +#define MIN_PER_HOUR 60
3.22 +#define SEC_PER_MIN 60
3.23 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
3.24 +
3.25 +static char *metadata_vars[] = {
3.26 + "lon",
3.27 + "lat",
3.28 + "lonrof",
3.29 + "latrof",
3.30 + "time",
3.31 + "hour", /* new metadata variable to be added to output files */
3.32 + "levsoi",
3.33 + "levlak",
3.34 + "edgen",
3.35 + "edgee",
3.36 + "edges",
3.37 + "edgew",
3.38 + "longxy",
3.39 + "latixy",
3.40 + "area",
3.41 + "areaupsc",
3.42 + "topo",
3.43 + "topodnsc",
3.44 + "landfrac",
3.45 + "landmask",
3.46 + "pftmask",
3.47 + "indxupsc",
3.48 + "mcdate",
3.49 + "mcsec",
3.50 + "mdcur",
3.51 + "mscur",
3.52 + "nstep",
3.53 + "time_bounds",
3.54 + "date_written",
3.55 + "time_written"
3.56 +};
3.57 +
3.58 +struct text_att {
3.59 + char *name;
3.60 + char *value;
3.61 + struct text_att *next;
3.62 + struct text_att *prev;
3.63 +};
3.64 +
3.65 +struct dim {
3.66 + int dimid;
3.67 + char *name;
3.68 + size_t len;
3.69 + struct dim *next;
3.70 + struct dim *prev;
3.71 +};
3.72 +
3.73 +struct var {
3.74 + int ncvarid;
3.75 + char *name;
3.76 + nc_type nctype;
3.77 + int ndims;
3.78 + int *dimids;
3.79 + int natts;
3.80 + char metadata;
3.81 + char FillFlag;
3.82 + float FillValue;
3.83 + struct var *next;
3.84 + struct var *prev;
3.85 +};
3.86 +
3.87 +static char *time_name = "time";
3.88 +static char *mcsec_name = "mcsec";
3.89 +static char *history_name = "history";
3.90 +static char *nsamples_name = "nsamples";
3.91 +static char *hour_name = "hour", *hour_long_name = "hour of day",
3.92 + *hour_units = "hour";
3.93 +static char *cell_method_name = "cell_method";
3.94 +/* input stuff */
3.95 +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
3.96 +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
3.97 +static struct text_att *input_gatt_head = NULL;
3.98 +static struct dim *input_dim_head = NULL, **input_dim_idx;
3.99 +static struct var *input_var_head = NULL;
3.100 +/* translation stuff */
3.101 +static int *idid2mdid, *idid2sdid; /* dimension IDs */
3.102 +/* output stuff */
3.103 +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
3.104 +static int mean_hour_dimid; /* special notes */
3.105 +static struct dim *mean_dim_head = NULL;
3.106 +static struct var *mean_var_head = NULL;
3.107 +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
3.108 +static int stddev_hour_dimid; /* special notes */
3.109 +static struct dim *stddev_dim_head = NULL;
3.110 +static struct var *stddev_var_head = NULL;
3.111 +
3.112 +char is_metadata(char *name)
3.113 +{
3.114 + int i;
3.115 +
3.116 + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
3.117 +
3.118 + if (i < nmvars)
3.119 + return 1;
3.120 + else
3.121 + return 0;
3.122 +}
3.123 +
3.124 +#if 0
3.125 +struct dim *dimlist_find_by_name(struct dim *head, char *name)
3.126 +{
3.127 + int i;
3.128 + struct dim *p = NULL;
3.129 +
3.130 + if (head) {
3.131 + node = head;
3.132 + for (i = 0 ; i == 0 || node != head; i++) {
3.133 + if (!strcmp(node->name, name))
3.134 + p = node;
3.135 + node = node->next;
3.136 + }
3.137 + return p;
3.138 + }
3.139 + else
3.140 + return NULL;
3.141 +}
3.142 +#endif
3.143 +
3.144 +struct var *varlist_find_by_name(struct var *head, char *name)
3.145 +{
3.146 + int i;
3.147 + struct var *node;
3.148 +
3.149 + if (head) {
3.150 + node = head;
3.151 + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
3.152 + node = node->next;
3.153 + if (i && node == head)
3.154 + return NULL;
3.155 + else
3.156 + return node;
3.157 + }
3.158 + else
3.159 + return NULL;
3.160 +}
3.161 +
3.162 +void gattlist_add_node(struct text_att **headp, char *name, char *value)
3.163 +{
3.164 + struct text_att *head, *node;
3.165 +
3.166 + head = *headp;
3.167 +
3.168 + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
3.169 + perror("gattlist_add_node");
3.170 + exit(2);
3.171 + }
3.172 +
3.173 + if (!(node->name = strdup(name))) {
3.174 + perror("gattlist_add_node");
3.175 + exit(2);
3.176 + }
3.177 + if (!(node->value = strdup(value))) {
3.178 + perror("gattlist_add_node");
3.179 + exit(2);
3.180 + }
3.181 +
3.182 + if (head == NULL) {
3.183 + /* first one */
3.184 + *headp = node;
3.185 + node->next = node;
3.186 + node->prev = node;
3.187 + }
3.188 + else {
3.189 + node->next = head;
3.190 + node->prev = head->prev;
3.191 + /* set this after setting node->prev from here */
3.192 + head->prev = node;
3.193 + /* set this after having set node->prev */
3.194 + node->prev->next = node;
3.195 + }
3.196 +
3.197 + return;
3.198 +}
3.199 +
3.200 +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
3.201 +{
3.202 + struct dim *head, *node;
3.203 +
3.204 + head = *headp;
3.205 +
3.206 + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
3.207 + perror("dimlist_add_node");
3.208 + exit(2);
3.209 + }
3.210 +
3.211 + node->dimid = dimid;
3.212 + if (!(node->name = strdup(name))) {
3.213 + perror("dimlist_add_node");
3.214 + exit(2);
3.215 + }
3.216 + node->len = len;
3.217 +
3.218 + if (head == NULL) {
3.219 + /* first one */
3.220 + *headp = node;
3.221 + node->next = node;
3.222 + node->prev = node;
3.223 + }
3.224 + else {
3.225 + node->next = head;
3.226 + node->prev = head->prev;
3.227 + /* set this after setting node->prev from here */
3.228 + head->prev = node;
3.229 + /* set this after having set node->prev */
3.230 + node->prev->next = node;
3.231 + }
3.232 +
3.233 + return node;
3.234 +}
3.235 +
3.236 +void varlist_add_node(struct var **headp, int ncvarid, char *name,
3.237 + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
3.238 + float FillValue)
3.239 +{
3.240 + int i;
3.241 + struct var *head, *node;
3.242 +
3.243 + head = *headp;
3.244 +
3.245 + if (!(node = (struct var *)malloc(sizeof(struct var)))) {
3.246 + perror("varlist_add_node");
3.247 + exit(3);
3.248 + }
3.249 +
3.250 + node->ncvarid = ncvarid;
3.251 + if (!(node->name = strdup(name))) {
3.252 + perror("varlist_add_node");
3.253 + exit(3);
3.254 + }
3.255 + node->nctype = nctype;
3.256 + node->ndims = ndims;
3.257 + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
3.258 + perror("varlist_add_node");
3.259 + exit(3);
3.260 + }
3.261 + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
3.262 + node->natts = natts;
3.263 + node->metadata = is_metadata(name);
3.264 + node->FillFlag = FillFlag;
3.265 + node->FillValue = FillValue;
3.266 +
3.267 + if (head == NULL) {
3.268 + /* first one */
3.269 + *headp = node;
3.270 + node->next = node;
3.271 + node->prev = node;
3.272 + }
3.273 + else {
3.274 + node->next = head;
3.275 + node->prev = head->prev;
3.276 + /* set this after setting node->prev from here */
3.277 + head->prev = node;
3.278 + /* set this after having set node->prev */
3.279 + node->prev->next = node;
3.280 + }
3.281 +
3.282 + return;
3.283 +}
3.284 +
3.285 +void varlist_print(struct var *headp)
3.286 +{
3.287 + int i, j;
3.288 + struct var *node;
3.289 +
3.290 + printf("Beginning of Variable List\n");
3.291 +
3.292 + if (headp) {
3.293 + node = headp;
3.294 + for (j = 0; j == 0 || node != headp; j++) {
3.295 + printf("Variable %d (ptr=%ld):\n", j, (long)node);
3.296 + printf("\tncvarid = %d\n", node->ncvarid);
3.297 + printf("\tname = %s\n", node->name);
3.298 + printf("\tnctype = %d\n", node->nctype);
3.299 + printf("\tndims = %d\n", node->ndims);
3.300 + printf("\tdimids =");
3.301 + for (i = 0; i < node->ndims; i++)
3.302 + printf(" %d", node->dimids[i]);
3.303 + printf("\n");
3.304 + printf("\tnatts = %d\n", node->natts);
3.305 + printf("\tmetadata = %d\n", (int)node->metadata);
3.306 + printf("\tnext = %ld\n", (long)node->next);
3.307 + printf("\tprev = %ld\n", (long)node->prev);
3.308 + node = node->next;
3.309 + }
3.310 + }
3.311 + else {
3.312 + printf("\tList undefined\n");
3.313 + }
3.314 +
3.315 + printf("End of Variable List\n");
3.316 +
3.317 + return;
3.318 +}
3.319 +
3.320 +void wrap_nc(int status)
3.321 +{
3.322 + if (status != NC_NOERR) {
3.323 + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
3.324 + exit(1);
3.325 + }
3.326 +
3.327 + return;
3.328 +}
3.329 +
3.330 +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
3.331 +{
3.332 + int i;
3.333 + char name[NC_MAX_NAME+1];
3.334 + size_t len;
3.335 + struct dim **in_dim_idx;
3.336 +
3.337 + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
3.338 + perror("*in_dim_idxp");
3.339 + exit(3);
3.340 + }
3.341 + in_dim_idx = *in_dim_idxp;
3.342 +
3.343 + for (i = 0; i < ndims; i++) {
3.344 + wrap_nc(nc_inq_dim(ncid, i, name, &len));
3.345 + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
3.346 + }
3.347 +
3.348 + return;
3.349 +}
3.350 +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
3.351 +{
3.352 + int i;
3.353 + char name[NC_MAX_NAME+1], *value;
3.354 + nc_type xtype;
3.355 + size_t len;
3.356 +
3.357 + for (i = 0; i < ngatts; i++) {
3.358 + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
3.359 + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
3.360 + if (xtype != NC_CHAR) {
3.361 + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
3.362 + exit(2);
3.363 + }
3.364 + if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
3.365 + perror("get_gatts: value");
3.366 + exit(3);
3.367 + }
3.368 + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
3.369 + value[(len+1-1)] = '\0';
3.370 + gattlist_add_node(gatt_headp, name, value);
3.371 + free(value); /* because gattlist_add_node() duplicates it */
3.372 + }
3.373 +
3.374 + return;
3.375 +}
3.376 +
3.377 +void get_vars(int ncid, int nvars, struct var **var_headp)
3.378 +{
3.379 + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
3.380 + size_t len;
3.381 + float FillValue;
3.382 + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
3.383 + nc_type xtype, atype;
3.384 +
3.385 + for (i = 0; i < nvars; i++) {
3.386 + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
3.387 + &natts));
3.388 + /* Check for _FillValue */
3.389 + FillFlag = 0;
3.390 + FillValue = 0.;
3.391 + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
3.392 + == NC_NOERR) {
3.393 + if (atype == NC_FLOAT && len) {
3.394 + wrap_nc(nc_get_att_float(ncid, i,
3.395 + fill_att_name, &FillValue));
3.396 + FillFlag = 1;
3.397 + }
3.398 + }
3.399 + varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
3.400 + natts, FillFlag, FillValue);
3.401 + }
3.402 +
3.403 + return;
3.404 +}
3.405 +
3.406 +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
3.407 + size_t nsamples, int **in2outp, int out_ncid,
3.408 + struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)
3.409 +{
3.410 + /*
3.411 + * Define dimensions on new files and build translation tables between
3.412 + * dimension IDs.
3.413 + */
3.414 + int j, dimid, ndims, *in2out;
3.415 + size_t len;
3.416 + struct dim *dnode;
3.417 +
3.418 + ndims = 0;
3.419 +
3.420 + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
3.421 + perror("put_dimensions: in2outp");
3.422 + exit(3);
3.423 + }
3.424 + in2out = *in2outp;
3.425 +
3.426 + if (in_dim_head) {
3.427 + dnode = in_dim_head;
3.428 + for (j = 0; j == 0 || dnode != in_dim_head; j++) {
3.429 + if (dnode->dimid == in_unlimdimid)
3.430 + len = NC_UNLIMITED;
3.431 + else
3.432 + len = dnode->len;
3.433 + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
3.434 + dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
3.435 + ++ndims;
3.436 + in2out[dnode->dimid] = dimid;
3.437 + if (dnode->dimid == in_unlimdimid)
3.438 + *out_unlimdimidp = dimid;
3.439 + /*
3.440 + * Just after the "time" dimension, add the new
3.441 + * "nsamples" and "hour" dimensions.
3.442 + */
3.443 + if (!strcmp(dnode->name, time_name)) {
3.444 + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
3.445 + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
3.446 + ++ndims;
3.447 +
3.448 + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
3.449 + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
3.450 + ++ndims;
3.451 + /* Track hour dimid for out files */
3.452 + *out_hour_dimidp = dimid;
3.453 + }
3.454 +
3.455 + dnode = dnode->next;
3.456 + }
3.457 + }
3.458 + else {
3.459 + fprintf(stderr, "WARNING: No dimensions defined!\n");
3.460 + }
3.461 +
3.462 + return ndims;
3.463 +}
3.464 +
3.465 +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
3.466 +{
3.467 + /*
3.468 + * Write out global attributes matching those of the input file.
3.469 + * Change history attribute to the string passed in as an argument.
3.470 + */
3.471 + int j, hflag = 0, ngatts;
3.472 + char *value;
3.473 + struct text_att *anode;
3.474 +
3.475 + ngatts = 0;
3.476 +
3.477 + if (in_gatt_head) {
3.478 + anode = in_gatt_head;
3.479 + for (j = 0; j == 0 || anode != in_gatt_head; j++) {
3.480 + if (!strcmp(anode->name, history_name)) {
3.481 + value = out_history;
3.482 + ++hflag;
3.483 + }
3.484 + else
3.485 + value = anode->value;
3.486 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
3.487 + ++ngatts;
3.488 + anode = anode->next;
3.489 + }
3.490 + /* If no history attribute on input, add one to the output */
3.491 + if (!hflag) {
3.492 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
3.493 + ++ngatts;
3.494 + }
3.495 + }
3.496 + else {
3.497 + fprintf(stderr, "WARNING: No global attributes defined!\n");
3.498 + }
3.499 +
3.500 + return ngatts;
3.501 +}
3.502 +
3.503 +int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)
3.504 +{
3.505 + /*
3.506 + * Translate between input dimension IDs and those for the output file.
3.507 + * For normal time-based variables, add a new dimension for hour of
3.508 + * day. For metadata variables, do not add new dimensions, even if
3.509 + * it is time-based. Return (possibly new) number of dimensions.
3.510 + */
3.511 + int i, ndims;
3.512 +
3.513 + for (i = ndims = 0; i < in_ndims; i++, ndims++) {
3.514 + out_dimids[ndims] = in2out[in_dimids[i]];
3.515 + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
3.516 + && !metadata) {
3.517 + ndims++;
3.518 + out_dimids[ndims] = hour_dimid;
3.519 + }
3.520 + }
3.521 +
3.522 + return ndims;
3.523 +}
3.524 +
3.525 +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
3.526 + int in_varid, int out_ncid, int out_varid)
3.527 +{
3.528 + /*
3.529 + * Copy the attributes of the input variable to those of the output
3.530 + * variable. If the variable is not a metadata variable, ensure that
3.531 + * the cell_method attribute is set either to "time: mean" or
3.532 + * "time: standard_deviation", depending on the output file type.
3.533 + */
3.534 +
3.535 + int i, natts;
3.536 + char name[NC_MAX_NAME+1], cmflag = 0;
3.537 + char *cell_methods[] = { "time: mean", "time: standard_deviation" };
3.538 +
3.539 + for (i = natts = 0; i < in_natts; i++, natts++) {
3.540 + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
3.541 + if (!strcmp(name, cell_method_name) && !metadata) {
3.542 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
3.543 + cmflag = 1;
3.544 + }
3.545 + else
3.546 + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
3.547 + }
3.548 + /*
3.549 + * If no cell_method attribute was specified for a non-metadata
3.550 + * variable on the input file, add the appropriate cell_method anyway
3.551 + * on the output file.
3.552 + */
3.553 + if (!cmflag && !metadata) {
3.554 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
3.555 + natts++;
3.556 + }
3.557 +
3.558 + return natts;
3.559 +}
3.560 +
3.561 +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
3.562 + int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp,
3.563 + int stat_type)
3.564 +{
3.565 + /*
3.566 + * Define variables on output file and copy attributes from input file.
3.567 + * Return number of variables defined.
3.568 + */
3.569 + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
3.570 + struct var *vnode;
3.571 +
3.572 + nvars = 0;
3.573 +
3.574 + if (in_var_head) {
3.575 + vnode = in_var_head;
3.576 + /*
3.577 + * March through input variables creating (mostly) the same
3.578 + * variables on the output file.
3.579 + */
3.580 + for (j = 0; j == 0 || vnode != in_var_head; j++) {
3.581 + ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
3.582 + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
3.583 + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
3.584 + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
3.585 + ++nvars;
3.586 + /*
3.587 + * Just after the "time" variable, add the new "hour"
3.588 + * variable for hour of day.
3.589 + */
3.590 + if (!strcmp(vnode->name, time_name)) {
3.591 + ndims = 1;
3.592 + dimids[0] = hour_dimid;
3.593 + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
3.594 + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
3.595 + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
3.596 + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
3.597 + ++nvars;
3.598 + }
3.599 +
3.600 + vnode = vnode->next;
3.601 + }
3.602 + }
3.603 + else {
3.604 + fprintf(stderr, "WARNING: No variables defined!\n");
3.605 + }
3.606 +
3.607 + return nvars;
3.608 +}
3.609 +
3.610 +void *alloc_var(nc_type nctype, size_t len)
3.611 +{
3.612 + void *val;
3.613 +
3.614 + switch (nctype) {
3.615 + case NC_FLOAT:
3.616 + if (!(val = (float *)malloc((len) * sizeof(float)))) {
3.617 + perror("alloc_var: val");
3.618 + exit(3);
3.619 + }
3.620 + break;
3.621 + case NC_INT:
3.622 + if (!(val = (int *)malloc((len) * sizeof(int)))) {
3.623 + perror("alloc_var: val");
3.624 + exit(3);
3.625 + }
3.626 + break;
3.627 + case NC_DOUBLE:
3.628 + if (!(val = (double *)malloc((len) * sizeof(double)))) {
3.629 + perror("alloc_var: val");
3.630 + exit(3);
3.631 + }
3.632 + break;
3.633 + case NC_CHAR:
3.634 + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
3.635 + perror("alloc_var: val");
3.636 + exit(3);
3.637 + }
3.638 + break;
3.639 + default:
3.640 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.641 + exit(3);
3.642 + }
3.643 +
3.644 + return val;
3.645 +}
3.646 +
3.647 +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
3.648 +{
3.649 + int i;
3.650 + size_t len = (size_t)1;
3.651 + void *val;
3.652 +
3.653 + /* Figure out the total size */
3.654 + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
3.655 +
3.656 + val = alloc_var(nctype, len);
3.657 +
3.658 + switch (nctype) {
3.659 + case NC_FLOAT:
3.660 + wrap_nc(nc_get_var_float(ncid, varid, val));
3.661 + break;
3.662 + case NC_INT:
3.663 + wrap_nc(nc_get_var_int(ncid, varid, val));
3.664 + break;
3.665 + case NC_DOUBLE:
3.666 + wrap_nc(nc_get_var_double(ncid, varid, val));
3.667 + break;
3.668 + case NC_CHAR:
3.669 + wrap_nc(nc_get_var_text(ncid, varid, val));
3.670 + break;
3.671 + default:
3.672 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.673 + exit(3);
3.674 + }
3.675 +
3.676 + return val;
3.677 +}
3.678 +
3.679 +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
3.680 +{
3.681 + int i;
3.682 + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
3.683 + void *val;
3.684 +
3.685 + /* Make sure time is really the first dimension */
3.686 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
3.687 + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
3.688 + exit(3);
3.689 + }
3.690 + /*
3.691 + * Figure out the total size for the slice (start at second dimension)
3.692 + * and build start/count vectors.
3.693 + */
3.694 + start[0] = tslice;
3.695 + count[0] = 1;
3.696 + for (i = 1; i < ndims; i++) {
3.697 + start[i] = 0;
3.698 + count[i] = (dim_idx[dimids[i]])->len;
3.699 + len *= (dim_idx[dimids[i]])->len;
3.700 + }
3.701 +
3.702 + val = alloc_var(nctype, len);
3.703 +
3.704 + switch (nctype) {
3.705 + case NC_FLOAT:
3.706 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
3.707 + break;
3.708 + case NC_INT:
3.709 + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
3.710 + break;
3.711 + case NC_DOUBLE:
3.712 + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
3.713 + break;
3.714 + case NC_CHAR:
3.715 + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
3.716 + break;
3.717 + default:
3.718 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.719 + exit(3);
3.720 + }
3.721 +
3.722 + return val;
3.723 +}
3.724 +
3.725 +void write_var(int ncid, int varid, nc_type nctype, void *val)
3.726 +{
3.727 + switch (nctype) {
3.728 + case NC_FLOAT:
3.729 + wrap_nc(nc_put_var_float(ncid, varid, val));
3.730 + break;
3.731 + case NC_INT:
3.732 + wrap_nc(nc_put_var_int(ncid, varid, val));
3.733 + break;
3.734 + case NC_DOUBLE:
3.735 + wrap_nc(nc_put_var_double(ncid, varid, val));
3.736 + break;
3.737 + case NC_CHAR:
3.738 + wrap_nc(nc_put_var_text(ncid, varid, val));
3.739 + break;
3.740 + default:
3.741 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.742 + exit(3);
3.743 + }
3.744 +
3.745 + return;
3.746 +}
3.747 +
3.748 +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
3.749 +{
3.750 + int i;
3.751 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
3.752 +
3.753 + /* Make sure time is really the first dimension */
3.754 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
3.755 + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
3.756 + exit(3);
3.757 + }
3.758 +
3.759 + /* Build start/count vectors */
3.760 + start[0] = tslice;
3.761 + count[0] = 1;
3.762 + for (i = 1; i < ndims; i++) {
3.763 + start[i] = 0;
3.764 + count[i] = (dim_idx[dimids[i]])->len;
3.765 + }
3.766 +
3.767 + switch (nctype) {
3.768 + case NC_FLOAT:
3.769 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
3.770 + break;
3.771 + case NC_INT:
3.772 + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
3.773 + break;
3.774 + case NC_DOUBLE:
3.775 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
3.776 + break;
3.777 + case NC_CHAR:
3.778 + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
3.779 + break;
3.780 + default:
3.781 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.782 + exit(3);
3.783 + }
3.784 +
3.785 + return;
3.786 +}
3.787 +
3.788 +void copy_metadata(int in_ncid, struct var *in_var_head,
3.789 + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head)
3.790 +{
3.791 + int i, j;
3.792 + float hr[HOURS_PER_DAY];
3.793 + struct var *in_vnode, *out_vnode;
3.794 + void *val;
3.795 +
3.796 + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
3.797 +
3.798 + if (in_var_head) {
3.799 + in_vnode = in_var_head;
3.800 + /*
3.801 + * March through input variables to stuff metadata values into
3.802 + * the new output files. NOTE: Time-based metadata variables
3.803 + * should contain only the last (time-slice) value (from all
3.804 + * input files).
3.805 + */
3.806 + for (j = 0; j == 0 || in_vnode != in_var_head; j++) {
3.807 + if (in_vnode->metadata) {
3.808 + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
3.809 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
3.810 + /* time is not the first dimension */
3.811 +#ifdef DEBUG
3.812 + printf("Copying metadata variable %s\n",
3.813 + in_vnode->name);
3.814 +#endif /* DEBUG */
3.815 + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
3.816 + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
3.817 + }
3.818 + else {
3.819 + /* time is the first dimension */
3.820 +#ifdef DEBUG
3.821 + printf("Copying last value of \
3.822 +time-based metadata variable %s\n", in_vnode->name);
3.823 +#endif /* DEBUG */
3.824 + val = read_timeslice(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, ((input_dim_idx[in_vnode->dimids[0]])->len - 1));
3.825 + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
3.826 + }
3.827 + free(val);
3.828 + /*
3.829 + * Just after the "time" variable, write out
3.830 + * the "hour" variable values.
3.831 + */
3.832 + if (!strcmp(in_vnode->name, time_name)) {
3.833 + out_vnode = varlist_find_by_name(out_var_head, hour_name);
3.834 + write_var(out_ncid, out_vnode->ncvarid,
3.835 + out_vnode->nctype, hr);
3.836 + }
3.837 + }
3.838 + else {
3.839 +#ifdef DEBUG
3.840 + printf("Skipping variable %s\n",
3.841 + in_vnode->name);
3.842 +#endif /* DEBUG */
3.843 + }
3.844 + in_vnode = in_vnode->next;
3.845 + }
3.846 + }
3.847 +
3.848 + return;
3.849 +}
3.850 +
3.851 +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
3.852 +{
3.853 + char *mean_history_gatt, *stddev_history_gatt,
3.854 + *mean_prefix="Statistical means from history files:",
3.855 + *stddev_prefix="Statistical standard deviations from history files:";
3.856 +
3.857 + /*
3.858 + * Construct strings for history global attributes for the two output
3.859 + * files.
3.860 + */
3.861 + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
3.862 + perror("open_inout:mean_history_gatt");
3.863 + exit(2);
3.864 + }
3.865 + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
3.866 + perror("open_inout:stddev_history_gatt");
3.867 + exit(2);
3.868 + }
3.869 + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
3.870 + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
3.871 +
3.872 + /* Open input file */
3.873 + wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid));
3.874 + /* Inquire about number of dimensions, variables, global attributes
3.875 + * and the ID of the unlimited dimension
3.876 + */
3.877 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
3.878 + &input_unlimdimid));
3.879 + /* Grab dimension IDs and lengths from input file */
3.880 + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
3.881 + /* Grab desired global attributes from input file */
3.882 + get_gatts(input_ncid, input_ngatts, &input_gatt_head);
3.883 + /* Grab variable IDs and attributes from input file */
3.884 + get_vars(input_ncid, input_nvars, &input_var_head);
3.885 +#ifdef DEBUG
3.886 + varlist_print(input_var_head);
3.887 +#endif /* DEBUG */
3.888 + /* Create netCDF files for monthly means and standard deviations */
3.889 + /* Create new netCDF files */
3.890 + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
3.891 + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
3.892 + /* Define dimensions */
3.893 + mean_ndims = put_dimensions(input_dim_head, input_ndims,
3.894 + input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
3.895 + &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
3.896 + stddev_ndims = put_dimensions(input_dim_head, input_ndims,
3.897 + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
3.898 + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
3.899 + /* Define variables and their attributes */
3.900 + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
3.901 + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
3.902 + MEAN_TYPE);
3.903 + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
3.904 + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
3.905 + STDDEV_TYPE);
3.906 + /* Store global attributes */
3.907 + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
3.908 + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
3.909 + stddev_history_gatt);
3.910 + /* End define mode */
3.911 + wrap_nc(nc_enddef(mean_ncid));
3.912 + wrap_nc(nc_enddef(stddev_ncid));
3.913 + /* Write out metdata variables that do not depend on "time" */
3.914 + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
3.915 + mean_var_head);
3.916 + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
3.917 + stddev_var_head);
3.918 +
3.919 + wrap_nc(nc_close(input_ncid));
3.920 +
3.921 + return;
3.922 +}
3.923 +
3.924 +size_t count_samples(int nifnames, char **input_fnames)
3.925 +{
3.926 + int i;
3.927 + char name[NC_MAX_NAME+1];
3.928 + size_t len, cnt = 0;
3.929 +
3.930 + for (i = 0; i < nifnames; i++) {
3.931 + /* Open input file */
3.932 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
3.933 + /*
3.934 + * Inquire about number of dimensions, variables, global
3.935 + * attributes and the ID of the unlimited dimension
3.936 + */
3.937 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
3.938 + &input_ngatts, &input_unlimdimid));
3.939 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
3.940 + if (strcmp(name, time_name)) {
3.941 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
3.942 + exit(4);
3.943 + }
3.944 +#ifdef DEBUG
3.945 + printf("%ld samples in %s\n", (long)len, input_fnames[i]);
3.946 +#endif /* DEBUG */
3.947 + wrap_nc(nc_close(input_ncid));
3.948 + cnt += len;
3.949 + }
3.950 + return cnt;
3.951 +}
3.952 +
3.953 +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
3.954 +{
3.955 + /*
3.956 + * Allocate space for arrays of means and standard deviations by
3.957 + * hour of day.
3.958 + */
3.959 + int i;
3.960 + size_t **cell_samples;
3.961 + double **means, **stddevs;
3.962 +
3.963 + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
3.964 + perror("alloc_means_stddevs:*meansp");
3.965 + exit(2);
3.966 + }
3.967 + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
3.968 + perror("alloc_means_stddevs:*stddevsp");
3.969 + exit(2);
3.970 + }
3.971 + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
3.972 + perror("alloc_means_stddevs:*cell_samplesp");
3.973 + exit(2);
3.974 + }
3.975 +
3.976 + means = *meansp;
3.977 + stddevs = *stddevsp;
3.978 + cell_samples = *cell_samplesp;
3.979 +
3.980 + for (i = 0; i < HOURS_PER_DAY; i++) {
3.981 + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
3.982 + perror("alloc_means_stddevs:means[i]");
3.983 + exit(2);
3.984 + }
3.985 + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
3.986 + perror("alloc_means_stddevs:stddevs[i]");
3.987 + exit(2);
3.988 + }
3.989 + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
3.990 + perror("alloc_means_stddevs:cell_samples[i]");
3.991 + exit(2);
3.992 + }
3.993 + }
3.994 +
3.995 + return;
3.996 +}
3.997 +
3.998 +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
3.999 +{
3.1000 + int hr, i, j, pos;
3.1001 +
3.1002 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
3.1003 + for (i = 0; i < d1; i++) {
3.1004 +#pragma _CRI concurrent
3.1005 + for (j = 0; j < d2; j++) {
3.1006 + pos = i*d2+j;
3.1007 + means[hr][pos] = FillValue;
3.1008 + stddevs[hr][pos] = FillValue;
3.1009 + cell_samples[hr][pos] = 0;
3.1010 + }
3.1011 + }
3.1012 + }
3.1013 +
3.1014 + return;
3.1015 +}
3.1016 +
3.1017 +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
3.1018 +{
3.1019 + int i, j, hr, pos;
3.1020 +
3.1021 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
3.1022 + for (i = 0; i < d1; i++) {
3.1023 +#pragma _CRI concurrent
3.1024 + for (j = 0; j < d2; j++) {
3.1025 + pos = i*d2+j;
3.1026 + cell_samples[hr][pos] = 0;
3.1027 + }
3.1028 + }
3.1029 + }
3.1030 +
3.1031 + return;
3.1032 +}
3.1033 +
3.1034 +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
3.1035 +{
3.1036 + /*
3.1037 + * Free space from arrays of means and standard deviations by
3.1038 + * hour of day.
3.1039 + */
3.1040 + int i;
3.1041 +
3.1042 + for (i = 0; i < HOURS_PER_DAY; i++) {
3.1043 + free(means[i]);
3.1044 + free(stddevs[i]);
3.1045 + free(cell_samples[i]);
3.1046 + }
3.1047 +
3.1048 + free(means);
3.1049 + free(stddevs);
3.1050 + free(cell_samples);
3.1051 +
3.1052 + return;
3.1053 +}
3.1054 +
3.1055 +void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
3.1056 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
3.1057 +{
3.1058 + int i, j, hr, pos;
3.1059 +
3.1060 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
3.1061 +
3.1062 + for (i = 0; i < d1; i++) {
3.1063 +#pragma _CRI concurrent
3.1064 + for (j = 0; j < d2; j++) {
3.1065 + pos = i*d2+j;
3.1066 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
3.1067 + if (cell_samples[hr][pos] == 0)
3.1068 + means[hr][pos] = (double)val[pos];
3.1069 + else
3.1070 + means[hr][pos] += (double)val[pos];
3.1071 + ++cell_samples[hr][pos];
3.1072 + }
3.1073 + }
3.1074 + }
3.1075 +
3.1076 + return;
3.1077 +}
3.1078 +
3.1079 +void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
3.1080 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
3.1081 +{
3.1082 + int i, j, hr, pos;
3.1083 +
3.1084 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
3.1085 +
3.1086 + for (i = 0; i < d1; i++) {
3.1087 +#pragma _CRI concurrent
3.1088 + for (j = 0; j < d2; j++) {
3.1089 + pos = i*d2+j;
3.1090 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
3.1091 + if (cell_samples[hr][pos] == 0)
3.1092 + means[hr][pos] = val[pos];
3.1093 + else
3.1094 + means[hr][pos] += val[pos];
3.1095 + ++cell_samples[hr][pos];
3.1096 + }
3.1097 + }
3.1098 + }
3.1099 +
3.1100 + return;
3.1101 +}
3.1102 +
3.1103 +
3.1104 +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
3.1105 +{
3.1106 + int i, j, hr, pos;
3.1107 +
3.1108 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
3.1109 + for (i = 0; i < d1; i++) {
3.1110 +#pragma _CRI concurrent
3.1111 + for (j = 0; j < d2; j++) {
3.1112 + pos = i*d2+j;
3.1113 + if (cell_samples[hr][pos] != 0) {
3.1114 + means[hr][pos] /= (double)cell_samples[hr][pos];
3.1115 + }
3.1116 + }
3.1117 + }
3.1118 + }
3.1119 +
3.1120 + return;
3.1121 +}
3.1122 +
3.1123 +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
3.1124 + char FillFlag, float FillValue, double **means, double **stddevs,
3.1125 + size_t **cell_samples)
3.1126 +{
3.1127 + int i, j, hr, pos;
3.1128 + double delta;
3.1129 +
3.1130 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
3.1131 +
3.1132 + for (i = 0; i < d1; i++) {
3.1133 +#pragma _CRI concurrent
3.1134 + for (j = 0; j < d2; j++) {
3.1135 + pos = i*d2+j;
3.1136 + if (!FillFlag || (FillFlag && val[pos] != FillValue
3.1137 + && means[hr][pos] != FillValue)) {
3.1138 + delta = means[hr][pos] - (double)val[pos];
3.1139 + if (cell_samples[hr][pos] == 0)
3.1140 + stddevs[hr][pos] = delta * delta;
3.1141 + else
3.1142 + stddevs[hr][pos] += delta * delta;
3.1143 + ++cell_samples[hr][pos];
3.1144 + }
3.1145 + }
3.1146 + }
3.1147 +
3.1148 + return;
3.1149 +}
3.1150 +
3.1151 +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
3.1152 + char FillFlag, float FillValue, double **means, double **stddevs,
3.1153 + size_t **cell_samples)
3.1154 +{
3.1155 + int i, j, hr, pos;
3.1156 + double delta;
3.1157 +
3.1158 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
3.1159 +
3.1160 + for (i = 0; i < d1; i++) {
3.1161 +#pragma _CRI concurrent
3.1162 + for (j = 0; j < d2; j++) {
3.1163 + pos = i*d2+j;
3.1164 + if (!FillFlag || (FillFlag && val[pos] != FillValue
3.1165 + && means[hr][pos] != FillValue)) {
3.1166 + delta = means[hr][pos] - val[pos];
3.1167 + if (cell_samples[hr][pos] == 0)
3.1168 + stddevs[hr][pos] = delta * delta;
3.1169 + else
3.1170 + stddevs[hr][pos] += delta * delta;
3.1171 + ++cell_samples[hr][pos];
3.1172 + }
3.1173 + }
3.1174 + }
3.1175 +
3.1176 + return;
3.1177 +}
3.1178 +
3.1179 +
3.1180 +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
3.1181 +{
3.1182 + int i, j, hr, pos;
3.1183 +
3.1184 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
3.1185 + for (i = 0; i < d1; i++) {
3.1186 +#pragma _CRI concurrent
3.1187 + for (j = 0; j < d2; j++) {
3.1188 + pos = i*d2+j;
3.1189 + if (cell_samples[hr][pos] != 0) {
3.1190 + stddevs[hr][pos] /= (double)cell_samples[hr][pos];
3.1191 + stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
3.1192 + }
3.1193 + }
3.1194 + }
3.1195 + }
3.1196 +
3.1197 + return;
3.1198 +}
3.1199 +
3.1200 +
3.1201 +void read_and_compute(int nifnames, char **input_fnames, size_t d1, size_t d2, char *var_name, nc_type nctype, int level, int ndims, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples)
3.1202 +{
3.1203 + int i, ts;
3.1204 + int varid, *mcsec;
3.1205 + size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
3.1206 + char name[NC_MAX_NAME+1];
3.1207 + void *val;
3.1208 +
3.1209 + /* Allocate space for timeslice */
3.1210 + val = alloc_var(nctype, (d1 * d2));
3.1211 +
3.1212 + for (i = 0; i < nifnames; i++) {
3.1213 +#ifdef DEBUG
3.1214 + printf("\tOpening %s", input_fnames[i]);
3.1215 + if (ndims > 3) printf(" to retrieve level %d", level);
3.1216 + if (stat_type == MEAN_TYPE)
3.1217 + printf(" for computing means\n");
3.1218 + else
3.1219 + printf(" for computing stddevs\n");
3.1220 +#endif /* DEBUG */
3.1221 + /* Open input file */
3.1222 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
3.1223 + /*
3.1224 + * Inquire about number of dimensions, variables, global
3.1225 + * attributes and the ID of the unlimited dimension
3.1226 + */
3.1227 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
3.1228 + &input_ngatts, &input_unlimdimid));
3.1229 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
3.1230 + if (strcmp(name, time_name)) {
3.1231 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
3.1232 + exit(4);
3.1233 + }
3.1234 + /* Get seconds of day variable */
3.1235 + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
3.1236 + if (!(mcsec = (int *)malloc(len * sizeof(int)))) {
3.1237 + perror("read_and_compute:mcsec");
3.1238 + exit(2);
3.1239 + }
3.1240 + wrap_nc(nc_get_var_int(input_ncid, varid, mcsec));
3.1241 + /* Get variable ID for requested variable */
3.1242 + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
3.1243 + /* Retrieve time slice of desired variable */
3.1244 + for (ts = 0; ts < len; ts++) {
3.1245 + start[0] = ts;
3.1246 + count[0] = 1;
3.1247 + start[(ndims-2)] = 0;
3.1248 + count[(ndims-2)] = d1;
3.1249 + start[(ndims-1)] = 0;
3.1250 + count[(ndims-1)] = d2;
3.1251 + if (ndims == 4) {
3.1252 + start[1] = level;
3.1253 + count[1] = 1;
3.1254 + }
3.1255 + switch(nctype) {
3.1256 + case NC_FLOAT:
3.1257 + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val));
3.1258 + if (stat_type == MEAN_TYPE)
3.1259 + add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
3.1260 + else
3.1261 + add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
3.1262 + break;
3.1263 + case NC_DOUBLE:
3.1264 + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val));
3.1265 + if (stat_type == MEAN_TYPE)
3.1266 + add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples);
3.1267 + else
3.1268 + add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
3.1269 + break;
3.1270 + default:
3.1271 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.1272 + exit(3);
3.1273 + }
3.1274 + }
3.1275 +
3.1276 + /* Free mcsec vector */
3.1277 + free(mcsec);
3.1278 + /* Close input file */
3.1279 + wrap_nc(nc_close(input_ncid));
3.1280 + }
3.1281 + /* Divide sums by number of samples */
3.1282 + if (stat_type == MEAN_TYPE)
3.1283 + divide_means(d1, d2, means, cell_samples);
3.1284 + else
3.1285 + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
3.1286 +
3.1287 + /* Free working variable array */
3.1288 + free(val);
3.1289 +
3.1290 + return;
3.1291 +}
3.1292 +
3.1293 +float *double_to_float(size_t len, double *dvar)
3.1294 +{
3.1295 + int i;
3.1296 + float *fvar;
3.1297 +
3.1298 + if (!(fvar = (float *)malloc(len * sizeof(float)))) {
3.1299 + perror("double_to_float:fvar");
3.1300 + exit(2);
3.1301 + }
3.1302 +
3.1303 + for (i = 0; i < len; i++)
3.1304 + fvar[i] = (float)dvar[i];
3.1305 +
3.1306 + return fvar;
3.1307 +}
3.1308 +
3.1309 +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
3.1310 + int level, int ndims, double **var_hours)
3.1311 +{
3.1312 + /* Output dimensions should be one larger than input dimensions */
3.1313 + int i, hr;
3.1314 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
3.1315 + float *fvar = NULL;
3.1316 +
3.1317 + if (nctype == NC_FLOAT) {
3.1318 + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
3.1319 + perror("write_var_hours:fvar");
3.1320 + exit(2);
3.1321 + }
3.1322 + }
3.1323 +
3.1324 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
3.1325 + start[0] = 0;
3.1326 + count[0] = 1; /* time */
3.1327 + start[1] = hr;
3.1328 + count[1] = 1; /* hour */
3.1329 + start[(ndims-2)] = 0;
3.1330 + count[(ndims-2)] = d1;
3.1331 + start[(ndims-1)] = 0;
3.1332 + count[(ndims-1)] = d2;
3.1333 + if (ndims == 5) {
3.1334 + start[2] = level;
3.1335 + count[2] = 1;
3.1336 + }
3.1337 + switch (nctype) {
3.1338 + case NC_FLOAT:
3.1339 + for (i = 0; i < (d1 * d2); i++)
3.1340 + fvar[i] = (float)var_hours[hr][i];
3.1341 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
3.1342 + break;
3.1343 + case NC_DOUBLE:
3.1344 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
3.1345 + break;
3.1346 + default:
3.1347 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
3.1348 + exit(3);
3.1349 + }
3.1350 + }
3.1351 +
3.1352 + if (nctype == NC_FLOAT)
3.1353 + free(fvar);
3.1354 +
3.1355 + return;
3.1356 +}
3.1357 +
3.1358 +void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
3.1359 +{
3.1360 + int j, k, nlevels;
3.1361 + size_t d1, d2, **cell_samples;
3.1362 + double **means, **stddevs;
3.1363 + struct var *in_vnode, *out_vnode;
3.1364 +
3.1365 + if (input_var_head) {
3.1366 + in_vnode = input_var_head;
3.1367 + /* March through non-metadata variables to compute stats */
3.1368 + for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
3.1369 + if (!in_vnode->metadata) {
3.1370 + /* Make sure time is really the first dimension */
3.1371 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
3.1372 + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
3.1373 + exit(3);
3.1374 + }
3.1375 + /* Figure out input dimensions */
3.1376 + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
3.1377 + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
3.1378 + exit(3);
3.1379 + }
3.1380 + /* Assume 2D output; find dimensions */
3.1381 + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
3.1382 + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
3.1383 + if (in_vnode->ndims == 3)
3.1384 + nlevels = 1;
3.1385 + else
3.1386 + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
3.1387 + /* Allocate working space for means and stddevs */
3.1388 + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
3.1389 + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
3.1390 + printf("Computing statistics for %s\n",
3.1391 + in_vnode->name);
3.1392 + /* Compute means and stddevs, then write them */
3.1393 + for (k = 0; k < nlevels; k++) {
3.1394 + /* Read and compute means */
3.1395 + read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
3.1396 + /* Find corresponding output variable on the mean output file */
3.1397 + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
3.1398 + /* Write out the means for this variable */
3.1399 + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
3.1400 + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
3.1401 + reset_cell_samples(d1, d2, cell_samples);
3.1402 + /* Read and compute stddevs using means */
3.1403 + read_and_compute(nifnames, input_fnames, d1, d2, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
3.1404 + /* Find corresponding output variable on the stddev output file */
3.1405 + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
3.1406 + /* Write out stddevs for this variable */
3.1407 + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
3.1408 + }
3.1409 +
3.1410 + /* Free working space for means and stddevs */
3.1411 + free_means_stddevs(means, stddevs, cell_samples);
3.1412 + }
3.1413 + in_vnode = in_vnode->next;
3.1414 + }
3.1415 + }
3.1416 + return;
3.1417 +}
3.1418 +
3.1419 +void usage(char *program)
3.1420 +{
3.1421 + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
3.1422 + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
3.1423 + exit(4);
3.1424 +}
3.1425 +
3.1426 +int main(int argc, char **argv)
3.1427 +{
3.1428 + int i, ic, nifnames;
3.1429 + size_t len, pos, nsamples;
3.1430 + char *mfname, *sfname, **ifnames, *flist;
3.1431 +
3.1432 + ifnames = NULL;
3.1433 + mfname = sfname = flist = NULL;
3.1434 + nifnames = 0;
3.1435 +
3.1436 +#ifdef DEBUG
3.1437 + printf("Number of metadata variables (nmvars) = %d\n", nmvars);
3.1438 + fflush(stdout);
3.1439 +#endif /* DEBUG */
3.1440 +
3.1441 + /* check command-line flags and switches */
3.1442 + while ((ic = getopt(argc, argv, "m:s:")) != -1) {
3.1443 + switch(ic) {
3.1444 + case 'm':
3.1445 + if (!(mfname = strdup(optarg))) {
3.1446 + perror("mfname");
3.1447 + exit(2);
3.1448 + }
3.1449 + break;
3.1450 + case 's':
3.1451 + if (!(sfname = strdup(optarg))) {
3.1452 + perror("sfname");
3.1453 + exit(2);
3.1454 + }
3.1455 + break;
3.1456 + }
3.1457 + }
3.1458 +
3.1459 + if (!mfname) {
3.1460 + fprintf(stderr, "Output file name for writing means required.\n");
3.1461 + usage(argv[0]);
3.1462 + }
3.1463 + if (!sfname) {
3.1464 + fprintf(stderr, "Output file name for writing standard deviations required.\n");
3.1465 + usage(argv[0]);
3.1466 + }
3.1467 +
3.1468 + if (optind < argc) {
3.1469 + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
3.1470 + perror("ifnames");
3.1471 + exit(2);
3.1472 + }
3.1473 + for (i = optind; i < argc; i++) {
3.1474 + if (!(ifnames[nifnames++] = strdup(argv[i]))) {
3.1475 + perror("ifnames[nifnames++]");
3.1476 + exit(2);
3.1477 + }
3.1478 + }
3.1479 + }
3.1480 + else {
3.1481 + fprintf(stderr, "At least one input file name is required.\n");
3.1482 + usage(argv[0]);
3.1483 + }
3.1484 +
3.1485 +
3.1486 + /*
3.1487 + * Build list of input files to be included in the global history
3.1488 + * attribute on the two outputs files.
3.1489 + */
3.1490 + for (i = len = 0; i < nifnames; i++)
3.1491 + len += strlen(ifnames[i]);
3.1492 + len += nifnames + 1; /* space in front of every name + null terminator */
3.1493 + if (!(flist = (char *)malloc(len * sizeof(char)))) {
3.1494 + perror("flist");
3.1495 + exit(2);
3.1496 + }
3.1497 + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
3.1498 + sprintf(&flist[pos], " %s", ifnames[i]);
3.1499 +#ifdef DEBUG
3.1500 + printf("flist=%s\n", flist);
3.1501 + fflush(stdout);
3.1502 +#endif /* DEBUG */
3.1503 +
3.1504 + /* Open every file to count up the number of time samples. */
3.1505 + nsamples = count_samples(nifnames, ifnames);
3.1506 +#ifdef DEBUG
3.1507 + printf("Number of samples across all files = %ld\n", (long)nsamples);
3.1508 + fflush(stdout);
3.1509 +#endif /* DEBUG */
3.1510 +
3.1511 + /*
3.1512 + * Open the *last* input file and the two output files (for mean and
3.1513 + * standard deviation). Define dimensions, variables, and attributes
3.1514 + * for the two output files based on the *last* input files. The
3.1515 + * last file is used because some metadata variables will contain
3.1516 + * only values from the last time slice from the period over which
3.1517 + * the statistics are computed.
3.1518 + */
3.1519 + open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples);
3.1520 +
3.1521 + compute_stats(nifnames, ifnames, nsamples);
3.1522 +
3.1523 + /* Close the two output files */
3.1524 + wrap_nc(nc_close(mean_ncid));
3.1525 + wrap_nc(nc_close(stddev_ncid));
3.1526 +
3.1527 + return 0;
3.1528 +}
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/h1_summary2.c Wed Sep 26 17:16:40 2007 -0400
4.3 @@ -0,0 +1,1614 @@
4.4 +#include <stdio.h>
4.5 +#include <unistd.h>
4.6 +#include <stdlib.h>
4.7 +#include <string.h>
4.8 +#include <math.h>
4.9 +#include "netcdf.h"
4.10 +
4.11 +/*
4.12 + * Loop through one month of hourly history tapes from the Community Land Model
4.13 + * and generate monthly statistics (means and standard deviations) of fields
4.14 + * by hour of day.
4.15 + */
4.16 +
4.17 +#define MEAN_TYPE 0
4.18 +#define STDDEV_TYPE 1
4.19 +
4.20 +#define HOURS_PER_DAY 24
4.21 +#define MIN_PER_HOUR 60
4.22 +#define SEC_PER_MIN 60
4.23 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN)
4.24 +
4.25 +static char *metadata_vars[] = {
4.26 + "lon",
4.27 + "lat",
4.28 + "lonrof",
4.29 + "latrof",
4.30 + "time",
4.31 + "hour", /* new metadata variable to be added to output files */
4.32 + "levsoi",
4.33 + "levlak",
4.34 + "edgen",
4.35 + "edgee",
4.36 + "edges",
4.37 + "edgew",
4.38 + "longxy",
4.39 + "latixy",
4.40 + "area",
4.41 + "areaupsc",
4.42 + "topo",
4.43 + "topodnsc",
4.44 + "landfrac",
4.45 + "landmask",
4.46 + "pftmask",
4.47 + "indxupsc",
4.48 + "mcdate",
4.49 + "mcsec",
4.50 + "mdcur",
4.51 + "mscur",
4.52 + "nstep",
4.53 + "time_bounds",
4.54 + "date_written",
4.55 + "time_written"
4.56 +};
4.57 +
4.58 +struct text_att {
4.59 + char *name;
4.60 + char *value;
4.61 + struct text_att *next;
4.62 + struct text_att *prev;
4.63 +};
4.64 +
4.65 +struct dim {
4.66 + int dimid;
4.67 + char *name;
4.68 + size_t len;
4.69 + struct dim *next;
4.70 + struct dim *prev;
4.71 +};
4.72 +
4.73 +struct var {
4.74 + int ncvarid;
4.75 + char *name;
4.76 + nc_type nctype;
4.77 + int ndims;
4.78 + int *dimids;
4.79 + int natts;
4.80 + char metadata;
4.81 + char FillFlag;
4.82 + float FillValue;
4.83 + struct var *next;
4.84 + struct var *prev;
4.85 +};
4.86 +
4.87 +static char *time_name = "time";
4.88 +static char *mcsec_name = "mcsec";
4.89 +static char *history_name = "history";
4.90 +static char *nsamples_name = "nsamples";
4.91 +static char *hour_name = "hour", *hour_long_name = "hour of day",
4.92 + *hour_units = "hour";
4.93 +static char *cell_method_name = "cell_method";
4.94 +/* input stuff */
4.95 +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars);
4.96 +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid;
4.97 +static struct text_att *input_gatt_head = NULL;
4.98 +static struct dim *input_dim_head = NULL, **input_dim_idx;
4.99 +static struct var *input_var_head = NULL;
4.100 +/* translation stuff */
4.101 +static int *idid2mdid, *idid2sdid; /* dimension IDs */
4.102 +/* output stuff */
4.103 +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid;
4.104 +static int mean_hour_dimid; /* special notes */
4.105 +static struct dim *mean_dim_head = NULL;
4.106 +static struct var *mean_var_head = NULL;
4.107 +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid;
4.108 +static int stddev_hour_dimid; /* special notes */
4.109 +static struct dim *stddev_dim_head = NULL;
4.110 +static struct var *stddev_var_head = NULL;
4.111 +
4.112 +char is_metadata(char *name)
4.113 +{
4.114 + int i;
4.115 +
4.116 + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++);
4.117 +
4.118 + if (i < nmvars)
4.119 + return 1;
4.120 + else
4.121 + return 0;
4.122 +}
4.123 +
4.124 +#if 0
4.125 +struct dim *dimlist_find_by_name(struct dim *head, char *name)
4.126 +{
4.127 + int i;
4.128 + struct dim *p = NULL;
4.129 +
4.130 + if (head) {
4.131 + node = head;
4.132 + for (i = 0 ; i == 0 || node != head; i++) {
4.133 + if (!strcmp(node->name, name))
4.134 + p = node;
4.135 + node = node->next;
4.136 + }
4.137 + return p;
4.138 + }
4.139 + else
4.140 + return NULL;
4.141 +}
4.142 +#endif
4.143 +
4.144 +struct var *varlist_find_by_name(struct var *head, char *name)
4.145 +{
4.146 + int i;
4.147 + struct var *node;
4.148 +
4.149 + if (head) {
4.150 + node = head;
4.151 + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++)
4.152 + node = node->next;
4.153 + if (i && node == head)
4.154 + return NULL;
4.155 + else
4.156 + return node;
4.157 + }
4.158 + else
4.159 + return NULL;
4.160 +}
4.161 +
4.162 +void gattlist_add_node(struct text_att **headp, char *name, char *value)
4.163 +{
4.164 + struct text_att *head, *node;
4.165 +
4.166 + head = *headp;
4.167 +
4.168 + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) {
4.169 + perror("gattlist_add_node");
4.170 + exit(2);
4.171 + }
4.172 +
4.173 + if (!(node->name = strdup(name))) {
4.174 + perror("gattlist_add_node");
4.175 + exit(2);
4.176 + }
4.177 + if (!(node->value = strdup(value))) {
4.178 + perror("gattlist_add_node");
4.179 + exit(2);
4.180 + }
4.181 +
4.182 + if (head == NULL) {
4.183 + /* first one */
4.184 + *headp = node;
4.185 + node->next = node;
4.186 + node->prev = node;
4.187 + }
4.188 + else {
4.189 + node->next = head;
4.190 + node->prev = head->prev;
4.191 + /* set this after setting node->prev from here */
4.192 + head->prev = node;
4.193 + /* set this after having set node->prev */
4.194 + node->prev->next = node;
4.195 + }
4.196 +
4.197 + return;
4.198 +}
4.199 +
4.200 +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len)
4.201 +{
4.202 + struct dim *head, *node;
4.203 +
4.204 + head = *headp;
4.205 +
4.206 + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) {
4.207 + perror("dimlist_add_node");
4.208 + exit(2);
4.209 + }
4.210 +
4.211 + node->dimid = dimid;
4.212 + if (!(node->name = strdup(name))) {
4.213 + perror("dimlist_add_node");
4.214 + exit(2);
4.215 + }
4.216 + node->len = len;
4.217 +
4.218 + if (head == NULL) {
4.219 + /* first one */
4.220 + *headp = node;
4.221 + node->next = node;
4.222 + node->prev = node;
4.223 + }
4.224 + else {
4.225 + node->next = head;
4.226 + node->prev = head->prev;
4.227 + /* set this after setting node->prev from here */
4.228 + head->prev = node;
4.229 + /* set this after having set node->prev */
4.230 + node->prev->next = node;
4.231 + }
4.232 +
4.233 + return node;
4.234 +}
4.235 +
4.236 +void varlist_add_node(struct var **headp, int ncvarid, char *name,
4.237 + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag,
4.238 + float FillValue)
4.239 +{
4.240 + int i;
4.241 + struct var *head, *node;
4.242 +
4.243 + head = *headp;
4.244 +
4.245 + if (!(node = (struct var *)malloc(sizeof(struct var)))) {
4.246 + perror("varlist_add_node");
4.247 + exit(3);
4.248 + }
4.249 +
4.250 + node->ncvarid = ncvarid;
4.251 + if (!(node->name = strdup(name))) {
4.252 + perror("varlist_add_node");
4.253 + exit(3);
4.254 + }
4.255 + node->nctype = nctype;
4.256 + node->ndims = ndims;
4.257 + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) {
4.258 + perror("varlist_add_node");
4.259 + exit(3);
4.260 + }
4.261 + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i];
4.262 + node->natts = natts;
4.263 + node->metadata = is_metadata(name);
4.264 + node->FillFlag = FillFlag;
4.265 + node->FillValue = FillValue;
4.266 +
4.267 + if (head == NULL) {
4.268 + /* first one */
4.269 + *headp = node;
4.270 + node->next = node;
4.271 + node->prev = node;
4.272 + }
4.273 + else {
4.274 + node->next = head;
4.275 + node->prev = head->prev;
4.276 + /* set this after setting node->prev from here */
4.277 + head->prev = node;
4.278 + /* set this after having set node->prev */
4.279 + node->prev->next = node;
4.280 + }
4.281 +
4.282 + return;
4.283 +}
4.284 +
4.285 +void varlist_print(struct var *headp)
4.286 +{
4.287 + int i, j;
4.288 + struct var *node;
4.289 +
4.290 + printf("Beginning of Variable List\n");
4.291 +
4.292 + if (headp) {
4.293 + node = headp;
4.294 + for (j = 0; j == 0 || node != headp; j++) {
4.295 + printf("Variable %d (ptr=%ld):\n", j, (long)node);
4.296 + printf("\tncvarid = %d\n", node->ncvarid);
4.297 + printf("\tname = %s\n", node->name);
4.298 + printf("\tnctype = %d\n", node->nctype);
4.299 + printf("\tndims = %d\n", node->ndims);
4.300 + printf("\tdimids =");
4.301 + for (i = 0; i < node->ndims; i++)
4.302 + printf(" %d", node->dimids[i]);
4.303 + printf("\n");
4.304 + printf("\tnatts = %d\n", node->natts);
4.305 + printf("\tmetadata = %d\n", (int)node->metadata);
4.306 + printf("\tnext = %ld\n", (long)node->next);
4.307 + printf("\tprev = %ld\n", (long)node->prev);
4.308 + node = node->next;
4.309 + }
4.310 + }
4.311 + else {
4.312 + printf("\tList undefined\n");
4.313 + }
4.314 +
4.315 + printf("End of Variable List\n");
4.316 +
4.317 + return;
4.318 +}
4.319 +
4.320 +void wrap_nc(int status)
4.321 +{
4.322 + if (status != NC_NOERR) {
4.323 + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
4.324 + exit(1);
4.325 + }
4.326 +
4.327 + return;
4.328 +}
4.329 +
4.330 +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp)
4.331 +{
4.332 + int i;
4.333 + char name[NC_MAX_NAME+1];
4.334 + size_t len;
4.335 + struct dim **in_dim_idx;
4.336 +
4.337 + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) {
4.338 + perror("*in_dim_idxp");
4.339 + exit(3);
4.340 + }
4.341 + in_dim_idx = *in_dim_idxp;
4.342 +
4.343 + for (i = 0; i < ndims; i++) {
4.344 + wrap_nc(nc_inq_dim(ncid, i, name, &len));
4.345 + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len);
4.346 + }
4.347 +
4.348 + return;
4.349 +}
4.350 +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp)
4.351 +{
4.352 + int i;
4.353 + char name[NC_MAX_NAME+1], *value;
4.354 + nc_type xtype;
4.355 + size_t len;
4.356 +
4.357 + for (i = 0; i < ngatts; i++) {
4.358 + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name));
4.359 + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len));
4.360 + if (xtype != NC_CHAR) {
4.361 + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name);
4.362 + exit(2);
4.363 + }
4.364 + if (!(value = (char *)malloc((len+1)*sizeof(char)))) {
4.365 + perror("get_gatts: value");
4.366 + exit(3);
4.367 + }
4.368 + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value));
4.369 + value[(len+1-1)] = '\0';
4.370 + gattlist_add_node(gatt_headp, name, value);
4.371 + free(value); /* because gattlist_add_node() duplicates it */
4.372 + }
4.373 +
4.374 + return;
4.375 +}
4.376 +
4.377 +void get_vars(int ncid, int nvars, struct var **var_headp)
4.378 +{
4.379 + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts;
4.380 + size_t len;
4.381 + float FillValue;
4.382 + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag;
4.383 + nc_type xtype, atype;
4.384 +
4.385 + for (i = 0; i < nvars; i++) {
4.386 + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids,
4.387 + &natts));
4.388 + /* Check for _FillValue */
4.389 + FillFlag = 0;
4.390 + FillValue = 0.;
4.391 + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len)
4.392 + == NC_NOERR) {
4.393 + if (atype == NC_FLOAT && len) {
4.394 + wrap_nc(nc_get_att_float(ncid, i,
4.395 + fill_att_name, &FillValue));
4.396 + FillFlag = 1;
4.397 + }
4.398 + }
4.399 + varlist_add_node(var_headp, i, name, xtype, ndims, dimids,
4.400 + natts, FillFlag, FillValue);
4.401 + }
4.402 +
4.403 + return;
4.404 +}
4.405 +
4.406 +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid,
4.407 + size_t nsamples, int **in2outp, int out_ncid,
4.408 + struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp)
4.409 +{
4.410 + /*
4.411 + * Define dimensions on new files and build translation tables between
4.412 + * dimension IDs.
4.413 + */
4.414 + int j, dimid, ndims, *in2out;
4.415 + size_t len;
4.416 + struct dim *dnode;
4.417 +
4.418 + ndims = 0;
4.419 +
4.420 + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) {
4.421 + perror("put_dimensions: in2outp");
4.422 + exit(3);
4.423 + }
4.424 + in2out = *in2outp;
4.425 +
4.426 + if (in_dim_head) {
4.427 + dnode = in_dim_head;
4.428 + for (j = 0; j == 0 || dnode != in_dim_head; j++) {
4.429 + if (dnode->dimid == in_unlimdimid)
4.430 + len = NC_UNLIMITED;
4.431 + else
4.432 + len = dnode->len;
4.433 + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid));
4.434 + dimlist_add_node(out_dim_headp, dimid, dnode->name, len);
4.435 + ++ndims;
4.436 + in2out[dnode->dimid] = dimid;
4.437 + if (dnode->dimid == in_unlimdimid)
4.438 + *out_unlimdimidp = dimid;
4.439 + /*
4.440 + * Just after the "time" dimension, add the new
4.441 + * "nsamples" and "hour" dimensions.
4.442 + */
4.443 + if (!strcmp(dnode->name, time_name)) {
4.444 + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid));
4.445 + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples);
4.446 + ++ndims;
4.447 +
4.448 + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid));
4.449 + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY);
4.450 + ++ndims;
4.451 + /* Track hour dimid for out files */
4.452 + *out_hour_dimidp = dimid;
4.453 + }
4.454 +
4.455 + dnode = dnode->next;
4.456 + }
4.457 + }
4.458 + else {
4.459 + fprintf(stderr, "WARNING: No dimensions defined!\n");
4.460 + }
4.461 +
4.462 + return ndims;
4.463 +}
4.464 +
4.465 +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history)
4.466 +{
4.467 + /*
4.468 + * Write out global attributes matching those of the input file.
4.469 + * Change history attribute to the string passed in as an argument.
4.470 + */
4.471 + int j, hflag = 0, ngatts;
4.472 + char *value;
4.473 + struct text_att *anode;
4.474 +
4.475 + ngatts = 0;
4.476 +
4.477 + if (in_gatt_head) {
4.478 + anode = in_gatt_head;
4.479 + for (j = 0; j == 0 || anode != in_gatt_head; j++) {
4.480 + if (!strcmp(anode->name, history_name)) {
4.481 + value = out_history;
4.482 + ++hflag;
4.483 + }
4.484 + else
4.485 + value = anode->value;
4.486 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value));
4.487 + ++ngatts;
4.488 + anode = anode->next;
4.489 + }
4.490 + /* If no history attribute on input, add one to the output */
4.491 + if (!hflag) {
4.492 + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history));
4.493 + ++ngatts;
4.494 + }
4.495 + }
4.496 + else {
4.497 + fprintf(stderr, "WARNING: No global attributes defined!\n");
4.498 + }
4.499 +
4.500 + return ngatts;
4.501 +}
4.502 +
4.503 +int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid)
4.504 +{
4.505 + /*
4.506 + * Translate between input dimension IDs and those for the output file.
4.507 + * For normal time-based variables, add a new dimension for hour of
4.508 + * day. For metadata variables, do not add new dimensions, even if
4.509 + * it is time-based. Return (possibly new) number of dimensions.
4.510 + */
4.511 + int i, ndims;
4.512 +
4.513 + for (i = ndims = 0; i < in_ndims; i++, ndims++) {
4.514 + out_dimids[ndims] = in2out[in_dimids[i]];
4.515 + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name)
4.516 + && !metadata) {
4.517 + ndims++;
4.518 + out_dimids[ndims] = hour_dimid;
4.519 + }
4.520 + }
4.521 +
4.522 + return ndims;
4.523 +}
4.524 +
4.525 +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid,
4.526 + int in_varid, int out_ncid, int out_varid)
4.527 +{
4.528 + /*
4.529 + * Copy the attributes of the input variable to those of the output
4.530 + * variable. If the variable is not a metadata variable, ensure that
4.531 + * the cell_method attribute is set either to "time: mean" or
4.532 + * "time: standard_deviation", depending on the output file type.
4.533 + */
4.534 +
4.535 + int i, natts;
4.536 + char name[NC_MAX_NAME+1], cmflag = 0;
4.537 + char *cell_methods[] = { "time: mean", "time: standard_deviation" };
4.538 +
4.539 + for (i = natts = 0; i < in_natts; i++, natts++) {
4.540 + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name));
4.541 + if (!strcmp(name, cell_method_name) && !metadata) {
4.542 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
4.543 + cmflag = 1;
4.544 + }
4.545 + else
4.546 + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid));
4.547 + }
4.548 + /*
4.549 + * If no cell_method attribute was specified for a non-metadata
4.550 + * variable on the input file, add the appropriate cell_method anyway
4.551 + * on the output file.
4.552 + */
4.553 + if (!cmflag && !metadata) {
4.554 + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type]));
4.555 + natts++;
4.556 + }
4.557 +
4.558 + return natts;
4.559 +}
4.560 +
4.561 +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head,
4.562 + int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp,
4.563 + int stat_type)
4.564 +{
4.565 + /*
4.566 + * Define variables on output file and copy attributes from input file.
4.567 + * Return number of variables defined.
4.568 + */
4.569 + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts;
4.570 + struct var *vnode;
4.571 +
4.572 + nvars = 0;
4.573 +
4.574 + if (in_var_head) {
4.575 + vnode = in_var_head;
4.576 + /*
4.577 + * March through input variables creating (mostly) the same
4.578 + * variables on the output file.
4.579 + */
4.580 + for (j = 0; j == 0 || vnode != in_var_head; j++) {
4.581 + ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid);
4.582 + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid));
4.583 + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid);
4.584 + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue);
4.585 + ++nvars;
4.586 + /*
4.587 + * Just after the "time" variable, add the new "hour"
4.588 + * variable for hour of day.
4.589 + */
4.590 + if (!strcmp(vnode->name, time_name)) {
4.591 + ndims = 1;
4.592 + dimids[0] = hour_dimid;
4.593 + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid));
4.594 + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name));
4.595 + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units));
4.596 + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0);
4.597 + ++nvars;
4.598 + }
4.599 +
4.600 + vnode = vnode->next;
4.601 + }
4.602 + }
4.603 + else {
4.604 + fprintf(stderr, "WARNING: No variables defined!\n");
4.605 + }
4.606 +
4.607 + return nvars;
4.608 +}
4.609 +
4.610 +void *alloc_var(nc_type nctype, size_t len)
4.611 +{
4.612 + void *val;
4.613 +
4.614 + switch (nctype) {
4.615 + case NC_FLOAT:
4.616 + if (!(val = (float *)malloc((len) * sizeof(float)))) {
4.617 + perror("alloc_var: val");
4.618 + exit(3);
4.619 + }
4.620 + break;
4.621 + case NC_INT:
4.622 + if (!(val = (int *)malloc((len) * sizeof(int)))) {
4.623 + perror("alloc_var: val");
4.624 + exit(3);
4.625 + }
4.626 + break;
4.627 + case NC_DOUBLE:
4.628 + if (!(val = (double *)malloc((len) * sizeof(double)))) {
4.629 + perror("alloc_var: val");
4.630 + exit(3);
4.631 + }
4.632 + break;
4.633 + case NC_CHAR:
4.634 + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) {
4.635 + perror("alloc_var: val");
4.636 + exit(3);
4.637 + }
4.638 + break;
4.639 + default:
4.640 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.641 + exit(3);
4.642 + }
4.643 +
4.644 + return val;
4.645 +}
4.646 +
4.647 +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx)
4.648 +{
4.649 + int i;
4.650 + size_t len = (size_t)1;
4.651 + void *val;
4.652 +
4.653 + /* Figure out the total size */
4.654 + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len;
4.655 +
4.656 + val = alloc_var(nctype, len);
4.657 +
4.658 + switch (nctype) {
4.659 + case NC_FLOAT:
4.660 + wrap_nc(nc_get_var_float(ncid, varid, val));
4.661 + break;
4.662 + case NC_INT:
4.663 + wrap_nc(nc_get_var_int(ncid, varid, val));
4.664 + break;
4.665 + case NC_DOUBLE:
4.666 + wrap_nc(nc_get_var_double(ncid, varid, val));
4.667 + break;
4.668 + case NC_CHAR:
4.669 + wrap_nc(nc_get_var_text(ncid, varid, val));
4.670 + break;
4.671 + default:
4.672 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.673 + exit(3);
4.674 + }
4.675 +
4.676 + return val;
4.677 +}
4.678 +
4.679 +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice)
4.680 +{
4.681 + int i;
4.682 + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
4.683 + void *val;
4.684 +
4.685 + /* Make sure time is really the first dimension */
4.686 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
4.687 + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name);
4.688 + exit(3);
4.689 + }
4.690 + /*
4.691 + * Figure out the total size for the slice (start at second dimension)
4.692 + * and build start/count vectors.
4.693 + */
4.694 + start[0] = tslice;
4.695 + count[0] = 1;
4.696 + for (i = 1; i < ndims; i++) {
4.697 + start[i] = 0;
4.698 + count[i] = (dim_idx[dimids[i]])->len;
4.699 + len *= (dim_idx[dimids[i]])->len;
4.700 + }
4.701 +
4.702 + val = alloc_var(nctype, len);
4.703 +
4.704 + switch (nctype) {
4.705 + case NC_FLOAT:
4.706 + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val));
4.707 + break;
4.708 + case NC_INT:
4.709 + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val));
4.710 + break;
4.711 + case NC_DOUBLE:
4.712 + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val));
4.713 + break;
4.714 + case NC_CHAR:
4.715 + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val));
4.716 + break;
4.717 + default:
4.718 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.719 + exit(3);
4.720 + }
4.721 +
4.722 + return val;
4.723 +}
4.724 +
4.725 +void write_var(int ncid, int varid, nc_type nctype, void *val)
4.726 +{
4.727 + switch (nctype) {
4.728 + case NC_FLOAT:
4.729 + wrap_nc(nc_put_var_float(ncid, varid, val));
4.730 + break;
4.731 + case NC_INT:
4.732 + wrap_nc(nc_put_var_int(ncid, varid, val));
4.733 + break;
4.734 + case NC_DOUBLE:
4.735 + wrap_nc(nc_put_var_double(ncid, varid, val));
4.736 + break;
4.737 + case NC_CHAR:
4.738 + wrap_nc(nc_put_var_text(ncid, varid, val));
4.739 + break;
4.740 + default:
4.741 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.742 + exit(3);
4.743 + }
4.744 +
4.745 + return;
4.746 +}
4.747 +
4.748 +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice)
4.749 +{
4.750 + int i;
4.751 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
4.752 +
4.753 + /* Make sure time is really the first dimension */
4.754 + if (strcmp((dim_idx[dimids[0]])->name, time_name)) {
4.755 + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name);
4.756 + exit(3);
4.757 + }
4.758 +
4.759 + /* Build start/count vectors */
4.760 + start[0] = tslice;
4.761 + count[0] = 1;
4.762 + for (i = 1; i < ndims; i++) {
4.763 + start[i] = 0;
4.764 + count[i] = (dim_idx[dimids[i]])->len;
4.765 + }
4.766 +
4.767 + switch (nctype) {
4.768 + case NC_FLOAT:
4.769 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val));
4.770 + break;
4.771 + case NC_INT:
4.772 + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val));
4.773 + break;
4.774 + case NC_DOUBLE:
4.775 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val));
4.776 + break;
4.777 + case NC_CHAR:
4.778 + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val));
4.779 + break;
4.780 + default:
4.781 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.782 + exit(3);
4.783 + }
4.784 +
4.785 + return;
4.786 +}
4.787 +
4.788 +void copy_metadata(int in_ncid, struct var *in_var_head,
4.789 + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head)
4.790 +{
4.791 + int i, j;
4.792 + float hr[HOURS_PER_DAY];
4.793 + struct var *in_vnode, *out_vnode;
4.794 + void *val;
4.795 +
4.796 + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i;
4.797 +
4.798 + if (in_var_head) {
4.799 + in_vnode = in_var_head;
4.800 + /*
4.801 + * March through input variables to stuff metadata values into
4.802 + * the new output files. NOTE: Time-based metadata variables
4.803 + * should contain only the last (time-slice) value (from all
4.804 + * input files).
4.805 + */
4.806 + for (j = 0; j == 0 || in_vnode != in_var_head; j++) {
4.807 + if (in_vnode->metadata) {
4.808 + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name);
4.809 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
4.810 + /* time is not the first dimension */
4.811 +#ifdef DEBUG
4.812 + printf("Copying metadata variable %s\n",
4.813 + in_vnode->name);
4.814 +#endif /* DEBUG */
4.815 + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx);
4.816 + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val);
4.817 + }
4.818 + else {
4.819 + /* time is the first dimension */
4.820 +#ifdef DEBUG
4.821 + printf("Copying last value of \
4.822 +time-based metadata variable %s\n", in_vnode->name);
4.823 +#endif /* DEBUG */
4.824 + val = read_timeslice(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, ((input_dim_idx[in_vnode->dimids[0]])->len - 1));
4.825 + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0);
4.826 + }
4.827 + free(val);
4.828 + /*
4.829 + * Just after the "time" variable, write out
4.830 + * the "hour" variable values.
4.831 + */
4.832 + if (!strcmp(in_vnode->name, time_name)) {
4.833 + out_vnode = varlist_find_by_name(out_var_head, hour_name);
4.834 + write_var(out_ncid, out_vnode->ncvarid,
4.835 + out_vnode->nctype, hr);
4.836 + }
4.837 + }
4.838 + else {
4.839 +#ifdef DEBUG
4.840 + printf("Skipping variable %s\n",
4.841 + in_vnode->name);
4.842 +#endif /* DEBUG */
4.843 + }
4.844 + in_vnode = in_vnode->next;
4.845 + }
4.846 + }
4.847 +
4.848 + return;
4.849 +}
4.850 +
4.851 +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples)
4.852 +{
4.853 + char *mean_history_gatt, *stddev_history_gatt,
4.854 + *mean_prefix="Statistical means from history files:",
4.855 + *stddev_prefix="Statistical standard deviations from history files:";
4.856 +
4.857 + /*
4.858 + * Construct strings for history global attributes for the two output
4.859 + * files.
4.860 + */
4.861 + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) {
4.862 + perror("open_inout:mean_history_gatt");
4.863 + exit(2);
4.864 + }
4.865 + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) {
4.866 + perror("open_inout:stddev_history_gatt");
4.867 + exit(2);
4.868 + }
4.869 + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist);
4.870 + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist);
4.871 +
4.872 + /* Open input file */
4.873 + wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid));
4.874 + /* Inquire about number of dimensions, variables, global attributes
4.875 + * and the ID of the unlimited dimension
4.876 + */
4.877 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts,
4.878 + &input_unlimdimid));
4.879 + /* Grab dimension IDs and lengths from input file */
4.880 + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx);
4.881 + /* Grab desired global attributes from input file */
4.882 + get_gatts(input_ncid, input_ngatts, &input_gatt_head);
4.883 + /* Grab variable IDs and attributes from input file */
4.884 + get_vars(input_ncid, input_nvars, &input_var_head);
4.885 +#ifdef DEBUG
4.886 + varlist_print(input_var_head);
4.887 +#endif /* DEBUG */
4.888 + /* Create netCDF files for monthly means and standard deviations */
4.889 + /* Create new netCDF files */
4.890 + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid));
4.891 + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid));
4.892 + /* Define dimensions */
4.893 + mean_ndims = put_dimensions(input_dim_head, input_ndims,
4.894 + input_unlimdimid, nsamples, &idid2mdid, mean_ncid,
4.895 + &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid);
4.896 + stddev_ndims = put_dimensions(input_dim_head, input_ndims,
4.897 + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid,
4.898 + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid);
4.899 + /* Define variables and their attributes */
4.900 + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
4.901 + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head,
4.902 + MEAN_TYPE);
4.903 + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head,
4.904 + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head,
4.905 + STDDEV_TYPE);
4.906 + /* Store global attributes */
4.907 + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt);
4.908 + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid,
4.909 + stddev_history_gatt);
4.910 + /* End define mode */
4.911 + wrap_nc(nc_enddef(mean_ncid));
4.912 + wrap_nc(nc_enddef(stddev_ncid));
4.913 + /* Write out metdata variables that do not depend on "time" */
4.914 + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid,
4.915 + mean_var_head);
4.916 + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid,
4.917 + stddev_var_head);
4.918 +
4.919 + wrap_nc(nc_close(input_ncid));
4.920 +
4.921 + return;
4.922 +}
4.923 +
4.924 +size_t count_samples(int nifnames, char **input_fnames)
4.925 +{
4.926 + int i;
4.927 + char name[NC_MAX_NAME+1];
4.928 + size_t len, cnt = 0;
4.929 +
4.930 + for (i = 0; i < nifnames; i++) {
4.931 + /* Open input file */
4.932 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
4.933 + /*
4.934 + * Inquire about number of dimensions, variables, global
4.935 + * attributes and the ID of the unlimited dimension
4.936 + */
4.937 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
4.938 + &input_ngatts, &input_unlimdimid));
4.939 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
4.940 + if (strcmp(name, time_name)) {
4.941 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
4.942 + exit(4);
4.943 + }
4.944 +#ifdef DEBUG
4.945 + printf("%ld samples in %s\n", (long)len, input_fnames[i]);
4.946 +#endif /* DEBUG */
4.947 + wrap_nc(nc_close(input_ncid));
4.948 + cnt += len;
4.949 + }
4.950 + return cnt;
4.951 +}
4.952 +
4.953 +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp)
4.954 +{
4.955 + /*
4.956 + * Allocate space for arrays of means and standard deviations by
4.957 + * hour of day.
4.958 + */
4.959 + int i;
4.960 + size_t **cell_samples;
4.961 + double **means, **stddevs;
4.962 +
4.963 + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
4.964 + perror("alloc_means_stddevs:*meansp");
4.965 + exit(2);
4.966 + }
4.967 + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) {
4.968 + perror("alloc_means_stddevs:*stddevsp");
4.969 + exit(2);
4.970 + }
4.971 + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) {
4.972 + perror("alloc_means_stddevs:*cell_samplesp");
4.973 + exit(2);
4.974 + }
4.975 +
4.976 + means = *meansp;
4.977 + stddevs = *stddevsp;
4.978 + cell_samples = *cell_samplesp;
4.979 +
4.980 + for (i = 0; i < HOURS_PER_DAY; i++) {
4.981 + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
4.982 + perror("alloc_means_stddevs:means[i]");
4.983 + exit(2);
4.984 + }
4.985 + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) {
4.986 + perror("alloc_means_stddevs:stddevs[i]");
4.987 + exit(2);
4.988 + }
4.989 + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) {
4.990 + perror("alloc_means_stddevs:cell_samples[i]");
4.991 + exit(2);
4.992 + }
4.993 + }
4.994 +
4.995 + return;
4.996 +}
4.997 +
4.998 +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue)
4.999 +{
4.1000 + int hr, i, j, pos;
4.1001 +
4.1002 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
4.1003 + for (i = 0; i < d1; i++) {
4.1004 +#pragma _CRI concurrent
4.1005 + for (j = 0; j < d2; j++) {
4.1006 + pos = i*d2+j;
4.1007 + means[hr][pos] = FillValue;
4.1008 + stddevs[hr][pos] = FillValue;
4.1009 + cell_samples[hr][pos] = 0;
4.1010 + }
4.1011 + }
4.1012 + }
4.1013 +
4.1014 + return;
4.1015 +}
4.1016 +
4.1017 +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples)
4.1018 +{
4.1019 + int i, j, hr, pos;
4.1020 +
4.1021 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
4.1022 + for (i = 0; i < d1; i++) {
4.1023 +#pragma _CRI concurrent
4.1024 + for (j = 0; j < d2; j++) {
4.1025 + pos = i*d2+j;
4.1026 + cell_samples[hr][pos] = 0;
4.1027 + }
4.1028 + }
4.1029 + }
4.1030 +
4.1031 + return;
4.1032 +}
4.1033 +
4.1034 +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples)
4.1035 +{
4.1036 + /*
4.1037 + * Free space from arrays of means and standard deviations by
4.1038 + * hour of day.
4.1039 + */
4.1040 + int i;
4.1041 +
4.1042 + for (i = 0; i < HOURS_PER_DAY; i++) {
4.1043 + free(means[i]);
4.1044 + free(stddevs[i]);
4.1045 + free(cell_samples[i]);
4.1046 + }
4.1047 +
4.1048 + free(means);
4.1049 + free(stddevs);
4.1050 + free(cell_samples);
4.1051 +
4.1052 + return;
4.1053 +}
4.1054 +
4.1055 +void add_to_means_float(float *val, int sec, size_t d1, size_t d2,
4.1056 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
4.1057 +{
4.1058 + int i, j, hr, pos;
4.1059 +
4.1060 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
4.1061 +
4.1062 + for (i = 0; i < d1; i++) {
4.1063 +#pragma _CRI concurrent
4.1064 + for (j = 0; j < d2; j++) {
4.1065 + pos = i*d2+j;
4.1066 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
4.1067 + if (cell_samples[hr][pos] == 0)
4.1068 + means[hr][pos] = (double)val[pos];
4.1069 + else
4.1070 + means[hr][pos] += (double)val[pos];
4.1071 + ++cell_samples[hr][pos];
4.1072 + }
4.1073 + }
4.1074 + }
4.1075 +
4.1076 + return;
4.1077 +}
4.1078 +
4.1079 +void add_to_means_double(double *val, int sec, size_t d1, size_t d2,
4.1080 + char FillFlag, float FillValue, double **means, size_t **cell_samples)
4.1081 +{
4.1082 + int i, j, hr, pos;
4.1083 +
4.1084 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
4.1085 +
4.1086 + for (i = 0; i < d1; i++) {
4.1087 +#pragma _CRI concurrent
4.1088 + for (j = 0; j < d2; j++) {
4.1089 + pos = i*d2+j;
4.1090 + if (!FillFlag || (FillFlag && val[pos] != FillValue)) {
4.1091 + if (cell_samples[hr][pos] == 0)
4.1092 + means[hr][pos] = val[pos];
4.1093 + else
4.1094 + means[hr][pos] += val[pos];
4.1095 + ++cell_samples[hr][pos];
4.1096 + }
4.1097 + }
4.1098 + }
4.1099 +
4.1100 + return;
4.1101 +}
4.1102 +
4.1103 +
4.1104 +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples)
4.1105 +{
4.1106 + int i, j, hr, pos;
4.1107 +
4.1108 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
4.1109 + for (i = 0; i < d1; i++) {
4.1110 +#pragma _CRI concurrent
4.1111 + for (j = 0; j < d2; j++) {
4.1112 + pos = i*d2+j;
4.1113 + if (cell_samples[hr][pos] != 0) {
4.1114 + means[hr][pos] /= (double)cell_samples[hr][pos];
4.1115 + }
4.1116 + }
4.1117 + }
4.1118 + }
4.1119 +
4.1120 + return;
4.1121 +}
4.1122 +
4.1123 +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2,
4.1124 + char FillFlag, float FillValue, double **means, double **stddevs,
4.1125 + size_t **cell_samples)
4.1126 +{
4.1127 + int i, j, hr, pos;
4.1128 + double delta;
4.1129 +
4.1130 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
4.1131 +
4.1132 + for (i = 0; i < d1; i++) {
4.1133 +#pragma _CRI concurrent
4.1134 + for (j = 0; j < d2; j++) {
4.1135 + pos = i*d2+j;
4.1136 + if (!FillFlag || (FillFlag && val[pos] != FillValue
4.1137 + && means[hr][pos] != FillValue)) {
4.1138 + delta = means[hr][pos] - (double)val[pos];
4.1139 + if (cell_samples[hr][pos] == 0)
4.1140 + stddevs[hr][pos] = delta * delta;
4.1141 + else
4.1142 + stddevs[hr][pos] += delta * delta;
4.1143 + ++cell_samples[hr][pos];
4.1144 + }
4.1145 + }
4.1146 + }
4.1147 +
4.1148 + return;
4.1149 +}
4.1150 +
4.1151 +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2,
4.1152 + char FillFlag, float FillValue, double **means, double **stddevs,
4.1153 + size_t **cell_samples)
4.1154 +{
4.1155 + int i, j, hr, pos;
4.1156 + double delta;
4.1157 +
4.1158 + hr = (int)floor((double)sec/(double)SEC_PER_HOUR);
4.1159 +
4.1160 + for (i = 0; i < d1; i++) {
4.1161 +#pragma _CRI concurrent
4.1162 + for (j = 0; j < d2; j++) {
4.1163 + pos = i*d2+j;
4.1164 + if (!FillFlag || (FillFlag && val[pos] != FillValue
4.1165 + && means[hr][pos] != FillValue)) {
4.1166 + delta = means[hr][pos] - val[pos];
4.1167 + if (cell_samples[hr][pos] == 0)
4.1168 + stddevs[hr][pos] = delta * delta;
4.1169 + else
4.1170 + stddevs[hr][pos] += delta * delta;
4.1171 + ++cell_samples[hr][pos];
4.1172 + }
4.1173 + }
4.1174 + }
4.1175 +
4.1176 + return;
4.1177 +}
4.1178 +
4.1179 +
4.1180 +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples)
4.1181 +{
4.1182 + int i, j, hr, pos;
4.1183 +
4.1184 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
4.1185 + for (i = 0; i < d1; i++) {
4.1186 +#pragma _CRI concurrent
4.1187 + for (j = 0; j < d2; j++) {
4.1188 + pos = i*d2+j;
4.1189 + if (cell_samples[hr][pos] != 0) {
4.1190 + stddevs[hr][pos] /= (double)cell_samples[hr][pos];
4.1191 + stddevs[hr][pos] = sqrt(stddevs[hr][pos]);
4.1192 + }
4.1193 + }
4.1194 + }
4.1195 + }
4.1196 +
4.1197 + return;
4.1198 +}
4.1199 +
4.1200 +void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp)
4.1201 +{
4.1202 + int i;
4.1203 + void **vals;
4.1204 +
4.1205 + /* Allocate space for all timeslices */
4.1206 + if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) {
4.1207 + perror("alloc_all_samples:*valsp");
4.1208 + exit(2);
4.1209 + }
4.1210 + vals = *valsp;
4.1211 + for (i = 0; i < nsamples; i++) {
4.1212 + vals[i] = alloc_var(nctype, (d1 * d2));
4.1213 + }
4.1214 + if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) {
4.1215 + perror("alloc_all_samples:*mcsecp");
4.1216 + exit(2);
4.1217 + }
4.1218 +
4.1219 + return;
4.1220 +}
4.1221 +
4.1222 +void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec)
4.1223 +{
4.1224 + int i, j;
4.1225 + float **fvals = NULL;
4.1226 + double **dvals = NULL;
4.1227 +
4.1228 + switch(nctype) {
4.1229 + case NC_FLOAT:
4.1230 + fvals = (float **)vals;
4.1231 + break;
4.1232 + case NC_DOUBLE:
4.1233 + dvals = (double **)vals;
4.1234 + break;
4.1235 + default:
4.1236 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.1237 + exit(3);
4.1238 + }
4.1239 +
4.1240 + for (i = 0; i < nsamples; i++) {
4.1241 + for (j = 0; j < (d1 * d2); j++) {
4.1242 + switch(nctype) {
4.1243 + case NC_FLOAT:
4.1244 + fvals[i][j] = (float)0.0;
4.1245 + break;
4.1246 + case NC_DOUBLE:
4.1247 + dvals[i][j] = (double)0.0;
4.1248 + break;
4.1249 + default:
4.1250 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.1251 + exit(3);
4.1252 + }
4.1253 + }
4.1254 + mcsec[i] = 0;
4.1255 + }
4.1256 +
4.1257 + return;
4.1258 +}
4.1259 +
4.1260 +void free_all_samples(size_t nsamples, void **vals, int *mcsec)
4.1261 +{
4.1262 + int i;
4.1263 +
4.1264 + for (i = 0; i < nsamples; i++)
4.1265 + if (vals[i]) free(vals[i]);
4.1266 +
4.1267 + free(vals);
4.1268 + free(mcsec);
4.1269 +
4.1270 + return;
4.1271 +}
4.1272 +
4.1273 +void read_all_samples(int nifnames, char **input_fnames, size_t d1, size_t d2, size_t nsamples, char *var_name, nc_type nctype, int level, int ndims, void **vals, int *mcsec)
4.1274 +{
4.1275 + int i, ts;
4.1276 + int varid;
4.1277 + size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
4.1278 + char name[NC_MAX_NAME+1];
4.1279 +
4.1280 + for (i = 0; i < nifnames; i++) {
4.1281 +#ifdef DEBUG
4.1282 + printf("\tOpening %s", input_fnames[i]);
4.1283 + if (ndims > 3) printf(" to retrieve level %d\n", level);
4.1284 +#endif /* DEBUG */
4.1285 + /* Open input file */
4.1286 + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid));
4.1287 + /*
4.1288 + * Inquire about number of dimensions, variables, global
4.1289 + * attributes and the ID of the unlimited dimension
4.1290 + */
4.1291 + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars,
4.1292 + &input_ngatts, &input_unlimdimid));
4.1293 + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len));
4.1294 + if (strcmp(name, time_name)) {
4.1295 + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]);
4.1296 + exit(4);
4.1297 + }
4.1298 + /* Make sure we don't run off the end of allocated space */
4.1299 + if ((nslice+len) > nsamples) {
4.1300 + fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples);
4.1301 + exit(4);
4.1302 + }
4.1303 + /* Get seconds of day variable */
4.1304 + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid));
4.1305 + wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice]));
4.1306 + /* Get variable ID for requested variable */
4.1307 + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid));
4.1308 + /* Retrieve time slice of desired variable */
4.1309 + for (ts = 0; ts < len; ts++) {
4.1310 + start[0] = ts;
4.1311 + count[0] = 1;
4.1312 + start[(ndims-2)] = 0;
4.1313 + count[(ndims-2)] = d1;
4.1314 + start[(ndims-1)] = 0;
4.1315 + count[(ndims-1)] = d2;
4.1316 + if (ndims == 4) {
4.1317 + start[1] = level;
4.1318 + count[1] = 1;
4.1319 + }
4.1320 + switch(nctype) {
4.1321 + case NC_FLOAT:
4.1322 + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice]));
4.1323 + break;
4.1324 + case NC_DOUBLE:
4.1325 + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice]));
4.1326 + break;
4.1327 + default:
4.1328 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.1329 + exit(3);
4.1330 + }
4.1331 + nslice++;
4.1332 + }
4.1333 +
4.1334 + /* Close input file */
4.1335 + wrap_nc(nc_close(input_ncid));
4.1336 + }
4.1337 +
4.1338 + return;
4.1339 +}
4.1340 +
4.1341 +void compute(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec, char FillFlag, float FillValue, int stat_type, double **means, double **stddevs, size_t **cell_samples)
4.1342 +{
4.1343 + int i;
4.1344 +
4.1345 + for (i = 0; i < nsamples; i++) {
4.1346 + switch(nctype) {
4.1347 + case NC_FLOAT:
4.1348 + if (stat_type == MEAN_TYPE)
4.1349 + add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
4.1350 + else
4.1351 + add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
4.1352 + break;
4.1353 + case NC_DOUBLE:
4.1354 + if (stat_type == MEAN_TYPE)
4.1355 + add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples);
4.1356 + else
4.1357 + add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples);
4.1358 + break;
4.1359 + default:
4.1360 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.1361 + exit(3);
4.1362 + }
4.1363 + }
4.1364 +
4.1365 + /* Divide sums by number of samples */
4.1366 + if (stat_type == MEAN_TYPE)
4.1367 + divide_means(d1, d2, means, cell_samples);
4.1368 + else
4.1369 + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples);
4.1370 +
4.1371 + return;
4.1372 +}
4.1373 +
4.1374 +float *double_to_float(size_t len, double *dvar)
4.1375 +{
4.1376 + int i;
4.1377 + float *fvar;
4.1378 +
4.1379 + if (!(fvar = (float *)malloc(len * sizeof(float)))) {
4.1380 + perror("double_to_float:fvar");
4.1381 + exit(2);
4.1382 + }
4.1383 +
4.1384 + for (i = 0; i < len; i++)
4.1385 + fvar[i] = (float)dvar[i];
4.1386 +
4.1387 + return fvar;
4.1388 +}
4.1389 +
4.1390 +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2,
4.1391 + int level, int ndims, double **var_hours)
4.1392 +{
4.1393 + /* Output dimensions should be one larger than input dimensions */
4.1394 + int i, hr;
4.1395 + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
4.1396 + float *fvar = NULL;
4.1397 +
4.1398 + if (nctype == NC_FLOAT) {
4.1399 + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) {
4.1400 + perror("write_var_hours:fvar");
4.1401 + exit(2);
4.1402 + }
4.1403 + }
4.1404 +
4.1405 + for (hr = 0; hr < HOURS_PER_DAY; hr++) {
4.1406 + start[0] = 0;
4.1407 + count[0] = 1; /* time */
4.1408 + start[1] = hr;
4.1409 + count[1] = 1; /* hour */
4.1410 + start[(ndims-2)] = 0;
4.1411 + count[(ndims-2)] = d1;
4.1412 + start[(ndims-1)] = 0;
4.1413 + count[(ndims-1)] = d2;
4.1414 + if (ndims == 5) {
4.1415 + start[2] = level;
4.1416 + count[2] = 1;
4.1417 + }
4.1418 + switch (nctype) {
4.1419 + case NC_FLOAT:
4.1420 + for (i = 0; i < (d1 * d2); i++)
4.1421 + fvar[i] = (float)var_hours[hr][i];
4.1422 + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar));
4.1423 + break;
4.1424 + case NC_DOUBLE:
4.1425 + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr]));
4.1426 + break;
4.1427 + default:
4.1428 + fprintf(stderr, "netCDF external data type %d not supported\n", nctype);
4.1429 + exit(3);
4.1430 + }
4.1431 + }
4.1432 +
4.1433 + if (nctype == NC_FLOAT)
4.1434 + free(fvar);
4.1435 +
4.1436 + return;
4.1437 +}
4.1438 +
4.1439 +void compute_stats(int nifnames, char **input_fnames, size_t nsamples)
4.1440 +{
4.1441 + int j, k, nlevels, *mcsec;
4.1442 + size_t d1, d2, **cell_samples;
4.1443 + double **means, **stddevs;
4.1444 + struct var *in_vnode, *out_vnode;
4.1445 + void **vals;
4.1446 +
4.1447 + if (input_var_head) {
4.1448 + in_vnode = input_var_head;
4.1449 + /* March through non-metadata variables to compute stats */
4.1450 + for (j = 0; j == 0 || in_vnode != input_var_head; j++) {
4.1451 + if (!in_vnode->metadata) {
4.1452 + /* Make sure time is really the first dimension */
4.1453 + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) {
4.1454 + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name);
4.1455 + exit(3);
4.1456 + }
4.1457 + /* Figure out input dimensions */
4.1458 + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) {
4.1459 + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims);
4.1460 + exit(3);
4.1461 + }
4.1462 + /* Assume 2D output; find dimensions */
4.1463 + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len;
4.1464 + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len;
4.1465 + if (in_vnode->ndims == 3)
4.1466 + nlevels = 1;
4.1467 + else
4.1468 + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len;
4.1469 + /* Allocate working space for means and stddevs */
4.1470 + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples);
4.1471 + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue);
4.1472 + /* Allocate and initize space for entire field across all files */
4.1473 + alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec);
4.1474 + init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec);
4.1475 + printf("Computing statistics for %s\n",
4.1476 + in_vnode->name);
4.1477 + /* Compute means and stddevs, then write them */
4.1478 + for (k = 0; k < nlevels; k++) {
4.1479 + /* Read the entire fields from all files */
4.1480 + read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec);
4.1481 + /* Compute means */
4.1482 + compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples);
4.1483 + /* Find corresponding output variable on the mean output file */
4.1484 + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name);
4.1485 + /* Write out the means for this variable */
4.1486 + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means);
4.1487 + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */
4.1488 + reset_cell_samples(d1, d2, cell_samples);
4.1489 + /* Compute stddevs using means */
4.1490 + compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples);
4.1491 + /* Find corresponding output variable on the stddev output file */
4.1492 + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name);
4.1493 + /* Write out stddevs for this variable */
4.1494 + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs);
4.1495 + }
4.1496 +
4.1497 + /* Free all samples */
4.1498 + free_all_samples(nsamples, vals, mcsec);
4.1499 + /* Free working space for means and stddevs */
4.1500 + free_means_stddevs(means, stddevs, cell_samples);
4.1501 + }
4.1502 + in_vnode = in_vnode->next;
4.1503 + }
4.1504 + }
4.1505 + return;
4.1506 +}
4.1507 +
4.1508 +void usage(char *program)
4.1509 +{
4.1510 + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program);
4.1511 + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n");
4.1512 + exit(4);
4.1513 +}
4.1514 +
4.1515 +int main(int argc, char **argv)
4.1516 +{
4.1517 + int i, ic, nifnames;
4.1518 + size_t len, pos, nsamples;
4.1519 + char *mfname, *sfname, **ifnames, *flist;
4.1520 +
4.1521 + ifnames = NULL;
4.1522 + mfname = sfname = flist = NULL;
4.1523 + nifnames = 0;
4.1524 +
4.1525 +#ifdef DEBUG
4.1526 + printf("Number of metadata variables (nmvars) = %d\n", nmvars);
4.1527 + fflush(stdout);
4.1528 +#endif /* DEBUG */
4.1529 +
4.1530 + /* check command-line flags and switches */
4.1531 + while ((ic = getopt(argc, argv, "m:s:")) != -1) {
4.1532 + switch(ic) {
4.1533 + case 'm':
4.1534 + if (!(mfname = strdup(optarg))) {
4.1535 + perror("mfname");
4.1536 + exit(2);
4.1537 + }
4.1538 + break;
4.1539 + case 's':
4.1540 + if (!(sfname = strdup(optarg))) {
4.1541 + perror("sfname");
4.1542 + exit(2);
4.1543 + }
4.1544 + break;
4.1545 + }
4.1546 + }
4.1547 +
4.1548 + if (!mfname) {
4.1549 + fprintf(stderr, "Output file name for writing means required.\n");
4.1550 + usage(argv[0]);
4.1551 + }
4.1552 + if (!sfname) {
4.1553 + fprintf(stderr, "Output file name for writing standard deviations required.\n");
4.1554 + usage(argv[0]);
4.1555 + }
4.1556 +
4.1557 + if (optind < argc) {
4.1558 + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) {
4.1559 + perror("ifnames");
4.1560 + exit(2);
4.1561 + }
4.1562 + for (i = optind; i < argc; i++) {
4.1563 + if (!(ifnames[nifnames++] = strdup(argv[i]))) {
4.1564 + perror("ifnames[nifnames++]");
4.1565 + exit(2);
4.1566 + }
4.1567 + }
4.1568 + }
4.1569 + else {
4.1570 + fprintf(stderr, "At least one input file name is required.\n");
4.1571 + usage(argv[0]);
4.1572 + }
4.1573 +
4.1574 +
4.1575 + /*
4.1576 + * Build list of input files to be included in the global history
4.1577 + * attribute on the two outputs files.
4.1578 + */
4.1579 + for (i = len = 0; i < nifnames; i++)
4.1580 + len += strlen(ifnames[i]);
4.1581 + len += nifnames + 1; /* space in front of every name + null terminator */
4.1582 + if (!(flist = (char *)malloc(len * sizeof(char)))) {
4.1583 + perror("flist");
4.1584 + exit(2);
4.1585 + }
4.1586 + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++)
4.1587 + sprintf(&flist[pos], " %s", ifnames[i]);
4.1588 +#ifdef DEBUG
4.1589 + printf("flist=%s\n", flist);
4.1590 + fflush(stdout);
4.1591 +#endif /* DEBUG */
4.1592 +
4.1593 + /* Open every file to count up the number of time samples. */
4.1594 + nsamples = count_samples(nifnames, ifnames);
4.1595 +#ifdef DEBUG
4.1596 + printf("Number of samples across all files = %ld\n", (long)nsamples);
4.1597 + fflush(stdout);
4.1598 +#endif /* DEBUG */
4.1599 +
4.1600 + /*
4.1601 + * Open the *last* input file and the two output files (for mean and
4.1602 + * standard deviation). Define dimensions, variables, and attributes
4.1603 + * for the two output files based on the *last* input files. The
4.1604 + * last file is used because some metadata variables will contain
4.1605 + * only values from the last time slice from the period over which
4.1606 + * the statistics are computed.
4.1607 + */
4.1608 + open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples);
4.1609 +
4.1610 + compute_stats(nifnames, ifnames, nsamples);
4.1611 +
4.1612 + /* Close the two output files */
4.1613 + wrap_nc(nc_close(mean_ncid));
4.1614 + wrap_nc(nc_close(stddev_ncid));
4.1615 +
4.1616 + return 0;
4.1617 +}
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/setup_phoenix.bash Wed Sep 26 17:16:40 2007 -0400
5.3 @@ -0,0 +1,17 @@
5.4 +#!/bin/bash
5.5 +source /opt/modules/modules/init/bash
5.6 +# for robin
5.7 +#module unload PrgEnv
5.8 +##module load pgi/6.1
5.9 +##module load netcdf_pgi/3.6.1
5.10 +#module load netcdf_gcc/3.6.1
5.11 +# for phoenix
5.12 +module purge
5.13 +module load open
5.14 +module load PrgEnv.5407
5.15 +module unload mpt
5.16 +module load mpt.2.4.0.6
5.17 +module load pbs
5.18 +module load netcdf/3.5.1_r4
5.19 +#
5.20 +module list
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/setup_robin1.bash Wed Sep 26 17:16:40 2007 -0400
6.3 @@ -0,0 +1,11 @@
6.4 +#!/bin/bash
6.5 +source /opt/modules/modules/init/bash
6.6 +# for robin1
6.7 +module unload PrgEnv
6.8 +#module load pgi/6.1
6.9 +#module load netcdf_pgi/3.6.1
6.10 +module load netcdf_gcc/3.6.1
6.11 +# for phoenix
6.12 +#module load netcdf/3.5.1_r4
6.13 +#
6.14 +module list