# HG changeset patch # User Forrest Hoffman # Date 1190841400 14400 # Node ID 3c02cce30be8fddfe7ca1a28dc959cc5bc330f8d Initial commit of post-processing codes for C-LAMP experiments add_total_fields - modifies model output netCDF files to add fields computed from fields within the files. h1_summary - computes means and standard deviations of hourly output netCDF files, creating two new netCDF files (one for the means and one for the standard deviations) for each month by hour of day. h1_summary2 - the same as h1_summary, but it uses more memory to read in more timeslices at once, so it may be faster on some machines. diff -r 000000000000 -r 3c02cce30be8 Makefile --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Makefile Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,38 @@ +CC=gcc +# +# robin1 and phoenix +LIBS=-L$(NETCDF)/lib -lnetcdf -lm +# Penguins +#LIBS=-L/usr/local/netcdf/netcdf-3.6.1-gcc+pgi/lib -lnetcdf -lm +# +# robin1 and phoenix +CPPFLAGS=-I$(NETCDF)/include +# Penguins +#CPPFLAGS=-I/usr/local/netcdf/netcdf-3.6.1-gcc+pgi/include +# phoenix +#CC=cc +#CFLAGS=-O -h list=m $(CPPFLAGS) +# robin1 and Penguins +CFLAGS=-g -Wall -O $(CPPFLAGS) + +all: h1_summary h1_summary2 add_total_fields + +h1_summary: h1_summary.o + $(CC) $(CFLAGS) -o $@ h1_summary.o $(LIBS) + +h1_summary2: h1_summary2.o + $(CC) $(CFLAGS) -o $@ h1_summary2.o $(LIBS) + +add_total_fields: add_total_fields.o + $(CC) $(CFLAGS) -o $@ add_total_fields.o $(LIBS) + +clean: + $(RM) -f h1_summary.o h1_summary + $(RM) -f h1_summary2.o h1_summary2 + $(RM) -f add_total_fields.o add_total_fields + +install: all + cp -p h1_summary $(HOME)/bin/ + cp -p h1_summary2 $(HOME)/bin/ + cp -p add_total_fields $(HOME)/bin/ + diff -r 000000000000 -r 3c02cce30be8 add_total_fields.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/add_total_fields.c Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,548 @@ +#include +#include +#include +#include +#include "netcdf.h" + +/* + * Loop through all history tapes from the Community Land Model adding new + * fields that are totals of the multi-level fields included in total_vars[] + * below. + */ + +static char *total_vars[] = { + "SOILICE", + "SOILLIQ" +}; +static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars); +static char *combo_vars[] = { + "LATENT", /* FCTR + FCEV + FGEV */ + "ET", /* QVEGE + QVEGT + QSOIL */ + "ALL_SKY_ALBEDO", /* FSR / FSDS */ + "BLACK_SKY_ALBEDO", /* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */ + "NETRAD", /* FSA - FIRA */ +}; +static char *combo_long_names[] = { + "latent heat flux", + "evapotranspiration", + "surface all-sky albedo", + "surface black-sky albedo", + "net radiation" +}; +static char *combo_units[] = { + "watt/m^2", + "m/s", + "unitless", + "unitless", + "watt/m^2" +}; +static float combo_FillValues[] = { + 1.e+36, + 1.e+36, + 1.e+36, + 1.e+36, + 1.e+36 +}; +static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars); +static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total "; +static char *time_name = "time", *lon_name = "lon", *lat_name = "lat"; +static char *long_name_name = "long_name", *units_name = "units", + *cell_method_name = "cell_method", *FillValue_name = "_FillValue", + *missing_value_name = "missing_value"; +static char *cell_method = "time: mean"; + +void wrap_nc(int status) +{ + if (status != NC_NOERR) { + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); + exit(1); + } + + return; +} + +void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1, + size_t d2, char FillFlag, float FillValue, size_t tslice, float *var, + double *tvar) +{ + int i, j; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + /* Initialize tvar to FillValue */ + for (i = 0; i < (d1 * d2); i++) + tvar[i] = (double)FillValue; + + + start[0] = tslice; + count[0] = 1; + start[2] = 0; + count[2] = d1; + start[3] = 0; + count[3] = d2; + + /* March through the levels, totalling as we go */ + for (j = 0; j < nlevels; j++) { + start[1] = j; + count[1] = 1; + wrap_nc(nc_get_vara_float(ncid, varid, start, count, var)); + for (i = 0; i < (d1 * d2); i++) { + if (tvar[i] == (double)FillValue) + tvar[i] = (double)var[i]; + else + tvar[i] += (double)var[i]; + } + } + + return; +} + +void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp) +{ + int i; + float **svars; + + if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) { + perror("alloc_svars:svarsp"); + exit(2); + } + svars = *svarsp; + for (i = 0; i < nsvars; i++) { + if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) { + perror("alloc_svars:svars[i]"); + exit(2); + } + } + return; +} + +void free_svars(size_t nsvars, float **svars) +{ + int i; + + for (i = 0; i < nsvars; i++) + free(svars[i]); + free(svars); + return; +} + +void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name, + float *var) +{ + int varid; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + start[0] = tslice; + count[0] = 1; + start[1] = 0; + count[1] = d1; + start[2] = 0; + count[2] = d2; + + /* Get variable id */ + wrap_nc(nc_inq_varid(ncid, name, &varid)); + /* Assume the correct size and read the variable */ + wrap_nc(nc_get_vara_float(ncid, varid, start, count, var)); + + return; +} + +void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, + int num, double *tvar) +{ + int i; + size_t nsvars; + float **svars; /* source variables */ + double denom; + + /* Initialize tvar to FillValue */ + for (i = 0; i < (d1 * d2); i++) + tvar[i] = (double)combo_FillValues[num]; + + switch (num) { + case 0: /* LATENT */ + nsvars = 3; + alloc_svars(nsvars, d1, d2, &svars); + /* read timeslices */ + read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]); + read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]); + read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]); + /* compute new variable */ + for (i = 0; i < (d1 * d2); i++) { + if (svars[0][i] != combo_FillValues[num] && + svars[1][i] != combo_FillValues[num] && + svars[2][i] != combo_FillValues[num]) { + tvar[i] = (double)svars[0][i] + + (double)svars[1][i] + + (double)svars[2][i]; + } + } + free_svars(nsvars, svars); + break; + case 1: /* ET */ + nsvars = 3; + alloc_svars(nsvars, d1, d2, &svars); + /* read timeslices */ + read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]); + read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]); + read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]); + /* compute new variable */ + for (i = 0; i < (d1 * d2); i++) { + if (svars[0][i] != combo_FillValues[num] && + svars[1][i] != combo_FillValues[num] && + svars[2][i] != combo_FillValues[num]) { + tvar[i] = (double)svars[0][i] + + (double)svars[1][i] + + (double)svars[2][i]; + } + } + free_svars(nsvars, svars); + break; + case 2: /* ALL_SKY_ALBEDO */ + nsvars = 2; + alloc_svars(nsvars, d1, d2, &svars); + /* read timeslices */ + read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]); + read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]); + /* compute new variable */ + for (i = 0; i < (d1 * d2); i++) { + if (svars[0][i] != combo_FillValues[num] && + svars[1][i] != combo_FillValues[num] && + (svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) { + tvar[i] = (double)svars[0][i] + / (double)svars[1][i]; + } + } + free_svars(nsvars, svars); + break; + case 3: /* BLACK_SKY_ALBEDO */ + nsvars = 4; + alloc_svars(nsvars, d1, d2, &svars); + /* read timeslices */ + read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]); + read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]); + read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]); + read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]); + /* compute new variable */ + for (i = 0; i < (d1 * d2); i++) { + if (svars[0][i] != combo_FillValues[num] && + svars[1][i] != combo_FillValues[num] && + svars[2][i] != combo_FillValues[num] && + svars[3][i] != combo_FillValues[num]) { + denom = (double)svars[2][i] + (double)svars[3][i]; + if (denom > 1.e-35 || denom < -1.e-35) { + tvar[i] = ((double)svars[0][i] + + (double)svars[1][i]) / + ((double)svars[2][i] + + (double)svars[3][i]); + } + } + } + free_svars(nsvars, svars); + break; + case 4: /* NETRAD */ + nsvars = 2; + alloc_svars(nsvars, d1, d2, &svars); + /* read timeslices */ + read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]); + read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]); + /* compute new variable */ + for (i = 0; i < (d1 * d2); i++) { + if (svars[0][i] != combo_FillValues[num] && + svars[1][i] != combo_FillValues[num]) { + tvar[i] = (double)svars[0][i] + - (double)svars[1][i]; + } + } + free_svars(nsvars, svars); + break; + default: + fprintf(stderr, "Unknown combination variable %d!\n", num); + exit(3); + } + + return; +} + +void write_timeslice(int ncid, int varid, size_t d1, size_t d2, + size_t tslice, float *var) +{ + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + start[0] = tslice; + count[0] = 1; + start[1] = 0; + count[1] = d1; + start[2] = 0; + count[2] = d2; + + wrap_nc(nc_put_vara_float(ncid, varid, start, count, var)); + + return; +} + +void usage(char *program) +{ + fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n", + program); + fprintf(stderr, "Requires at least one history output file name\n"); + exit(5); +} + +int main(int argc, char **argv) +{ + char **ifnames; + int i, j, k, nifnames; + int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts, + dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS], + new_varid, lat_dimid, lon_dimid; + nc_type xtype, atype; + char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1], + att_name[NC_MAX_NAME+1], *longname, *new_longname; + size_t tlen, ts, alen, nlevels, d1, d2; + float FillValue; + char FillFlag; + float *var; + double *tvar; + +#ifdef DEBUG + printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars); + printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars); +#endif /* DEBUG */ + + /* Require at least one parameter (the file on which to work) */ + if (argc < 2) + usage(argv[0]); + + /* Put the list of filenames in an array of strings */ + if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) { + perror("ifnames"); + exit(2); + } + nifnames = 0; + for (i = 1; i < argc; i++) { + if (!(ifnames[nifnames++] = strdup(argv[i]))) { + perror("ifnames[nifnames++]"); + exit(2); + } + } + + /* + * March through all input files adding new variables for totals + * across levels. + */ + for (i = 0; i < nifnames; i++) { + printf("Processing %s\n", ifnames[i]); + /* Open file */ + wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid)); + wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]); + exit(4); + } + /* Retrieve lat and lon dimension IDs */ + wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid)); + wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid)); + /* Put file into define mode */ + wrap_nc(nc_redef(ncid)); + /* + * Define all of the new variables first to improve + * performance. This means some calls are made twice. + */ + /* First, add the new total variables */ + for (j = 0; j < nmtotal_vars; j++) { + /* Figure out source variable information */ + wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid)); + wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts)); + /* Make sure it has the correct number of dimensions */ + if (ndims != 4) { + fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims); + exit(4); + } + /* Make sure time is the first dimension */ + if (dimids[0] != unlimdimid) { + fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name); + exit(4); + } + /* Make sure it is a float type */ + if (xtype != NC_FLOAT) { + fprintf(stderr, "Only NC_FLOAT type is supported presently!\n"); + exit(4); + } + /* Set dimensions for new variable */ + new_ndims = 3; + new_dimids[0] = unlimdimid; + new_dimids[1] = dimids[2]; + new_dimids[2] = dimids[3]; + /* Define new variable */ + sprintf(new_name, "%s%s", varname_total_prefix, name); + printf("\tAdding variable %s\n", new_name); + wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims, + new_dimids, &new_varid)); + /* Copy attributes from source variable to new one */ + for (k = 0; k < natts; k++) { + wrap_nc(nc_inq_attname(ncid, varid, k, att_name)); + if (!strcmp(att_name, "long_name")) { + wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen)); + if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) { + perror("longname"); + exit(2); + } + wrap_nc(nc_get_att_text(ncid, varid, att_name, longname)); + longname[alen] = '\0'; + if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) { + perror("new_longname"); + exit(2); + } + sprintf(new_longname, "%s%s", longname_total_prefix, longname); + wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname)); + free(new_longname); + free(longname); + } + else + wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid)); + } + } + /* Second, add the new total variables */ + for (j = 0; j < nmcombo_vars; j++) { + /* Set dimensions for new variable */ + new_ndims = 3; + new_dimids[0] = unlimdimid; + new_dimids[1] = lat_dimid; + new_dimids[2] = lon_dimid; + /* Define new variable */ + printf("\tAdding variable %s\n", combo_vars[j]); + wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT, + new_ndims, new_dimids, &new_varid)); + /* Add attributes */ + wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name, + strlen(combo_long_names[j]), + combo_long_names[j])); + wrap_nc(nc_put_att_text(ncid, new_varid, units_name, + strlen(combo_units[j]), combo_units[j])); + wrap_nc(nc_put_att_text(ncid, new_varid, + cell_method_name, strlen(cell_method), + cell_method)); + wrap_nc(nc_put_att_float(ncid, new_varid, + FillValue_name, NC_FLOAT, 1, + &combo_FillValues[j])); + wrap_nc(nc_put_att_float(ncid, new_varid, + missing_value_name, NC_FLOAT, 1, + &combo_FillValues[j])); + } + /* Leave define mode */ + wrap_nc(nc_enddef(ncid)); + /* Now actually compute and write out the new total variables */ + for (j = 0; j < nmtotal_vars; j++) { + /* Figure out source variable information */ + wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid)); + wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts)); + /* Check for _FillValue */ + FillFlag = 0; + FillValue = 0.; + if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) { + if (atype == NC_FLOAT && alen) { + wrap_nc(nc_get_att_float(ncid, varid, + FillValue_name, &FillValue)); + FillFlag = 1; + } + } + /* Set dimensions for new variable */ + new_ndims = 3; + new_dimids[0] = unlimdimid; + new_dimids[1] = dimids[2]; + new_dimids[2] = dimids[3]; + + sprintf(new_name, "%s%s", varname_total_prefix, name); + printf("\tComputing and writing %s\n", new_name); + /* Retrieve the new varid again */ + wrap_nc(nc_inq_varid(ncid, new_name, &new_varid)); + /* + * Figure out dimensions and total size + * for a single time slice + */ + wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels)); + wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1)); + wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2)); + + /* + * Allocate space for reading in time slice and for + * the sum + */ + if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) { + perror("var"); + exit(2); + } + if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) { + perror("tvar"); + exit(2); + } + + /* Read timeslices, total, and write out new timeslices */ + for (ts = 0; ts < tlen; ts++) { + read_total_timeslice(ncid, varid, + nlevels, d1, d2, FillFlag, FillValue, + ts, var, tvar); + for (k = 0; k < (d1*d2); k++) + var[k] = (float)tvar[k]; + write_timeslice(ncid, new_varid, d1, d2, + ts, var); + } + + free(var); + free(tvar); + } + /* Now actually compute and write out the new combo variables */ + for (j = 0; j < nmcombo_vars; j++) { + /* Set dimensions for new variable */ + new_ndims = 3; + new_dimids[0] = unlimdimid; + new_dimids[1] = lat_dimid; + new_dimids[2] = lon_dimid; + + printf("\tComputing and writing %s\n", combo_vars[j]); + /* Retrieve the new varid again */ + wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid)); + /* + * Figure out dimensions and total size + * for a single time slice + */ + wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1)); + wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2)); + + /* + * Allocate space for reading in time slice and for + * the sum + */ + if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) { + perror("var"); + exit(2); + } + if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) { + perror("tvar"); + exit(2); + } + + /* Read timeslices, compute new variable, and write */ + for (ts = 0; ts < tlen; ts++) { + read_compute_timeslice(ncid, d1, d2, ts, + j, tvar); + for (k = 0; k < (d1*d2); k++) + var[k] = (float)tvar[k]; + write_timeslice(ncid, new_varid, d1, d2, + ts, var); + } + + free(var); + free(tvar); + } + /* Close file */ + wrap_nc(nc_close(ncid)); + } + + return 0; +} diff -r 000000000000 -r 3c02cce30be8 h1_summary.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/h1_summary.c Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,1525 @@ +#include +#include +#include +#include +#include +#include "netcdf.h" + +/* + * Loop through one month of hourly history tapes from the Community Land Model + * and generate monthly statistics (means and standard deviations) of fields + * by hour of day. + */ + +#define MEAN_TYPE 0 +#define STDDEV_TYPE 1 + +#define HOURS_PER_DAY 24 +#define MIN_PER_HOUR 60 +#define SEC_PER_MIN 60 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN) + +static char *metadata_vars[] = { + "lon", + "lat", + "lonrof", + "latrof", + "time", + "hour", /* new metadata variable to be added to output files */ + "levsoi", + "levlak", + "edgen", + "edgee", + "edges", + "edgew", + "longxy", + "latixy", + "area", + "areaupsc", + "topo", + "topodnsc", + "landfrac", + "landmask", + "pftmask", + "indxupsc", + "mcdate", + "mcsec", + "mdcur", + "mscur", + "nstep", + "time_bounds", + "date_written", + "time_written" +}; + +struct text_att { + char *name; + char *value; + struct text_att *next; + struct text_att *prev; +}; + +struct dim { + int dimid; + char *name; + size_t len; + struct dim *next; + struct dim *prev; +}; + +struct var { + int ncvarid; + char *name; + nc_type nctype; + int ndims; + int *dimids; + int natts; + char metadata; + char FillFlag; + float FillValue; + struct var *next; + struct var *prev; +}; + +static char *time_name = "time"; +static char *mcsec_name = "mcsec"; +static char *history_name = "history"; +static char *nsamples_name = "nsamples"; +static char *hour_name = "hour", *hour_long_name = "hour of day", + *hour_units = "hour"; +static char *cell_method_name = "cell_method"; +/* input stuff */ +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars); +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid; +static struct text_att *input_gatt_head = NULL; +static struct dim *input_dim_head = NULL, **input_dim_idx; +static struct var *input_var_head = NULL; +/* translation stuff */ +static int *idid2mdid, *idid2sdid; /* dimension IDs */ +/* output stuff */ +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid; +static int mean_hour_dimid; /* special notes */ +static struct dim *mean_dim_head = NULL; +static struct var *mean_var_head = NULL; +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid; +static int stddev_hour_dimid; /* special notes */ +static struct dim *stddev_dim_head = NULL; +static struct var *stddev_var_head = NULL; + +char is_metadata(char *name) +{ + int i; + + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++); + + if (i < nmvars) + return 1; + else + return 0; +} + +#if 0 +struct dim *dimlist_find_by_name(struct dim *head, char *name) +{ + int i; + struct dim *p = NULL; + + if (head) { + node = head; + for (i = 0 ; i == 0 || node != head; i++) { + if (!strcmp(node->name, name)) + p = node; + node = node->next; + } + return p; + } + else + return NULL; +} +#endif + +struct var *varlist_find_by_name(struct var *head, char *name) +{ + int i; + struct var *node; + + if (head) { + node = head; + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++) + node = node->next; + if (i && node == head) + return NULL; + else + return node; + } + else + return NULL; +} + +void gattlist_add_node(struct text_att **headp, char *name, char *value) +{ + struct text_att *head, *node; + + head = *headp; + + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) { + perror("gattlist_add_node"); + exit(2); + } + + if (!(node->name = strdup(name))) { + perror("gattlist_add_node"); + exit(2); + } + if (!(node->value = strdup(value))) { + perror("gattlist_add_node"); + exit(2); + } + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len) +{ + struct dim *head, *node; + + head = *headp; + + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) { + perror("dimlist_add_node"); + exit(2); + } + + node->dimid = dimid; + if (!(node->name = strdup(name))) { + perror("dimlist_add_node"); + exit(2); + } + node->len = len; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return node; +} + +void varlist_add_node(struct var **headp, int ncvarid, char *name, + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag, + float FillValue) +{ + int i; + struct var *head, *node; + + head = *headp; + + if (!(node = (struct var *)malloc(sizeof(struct var)))) { + perror("varlist_add_node"); + exit(3); + } + + node->ncvarid = ncvarid; + if (!(node->name = strdup(name))) { + perror("varlist_add_node"); + exit(3); + } + node->nctype = nctype; + node->ndims = ndims; + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) { + perror("varlist_add_node"); + exit(3); + } + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i]; + node->natts = natts; + node->metadata = is_metadata(name); + node->FillFlag = FillFlag; + node->FillValue = FillValue; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +void varlist_print(struct var *headp) +{ + int i, j; + struct var *node; + + printf("Beginning of Variable List\n"); + + if (headp) { + node = headp; + for (j = 0; j == 0 || node != headp; j++) { + printf("Variable %d (ptr=%ld):\n", j, (long)node); + printf("\tncvarid = %d\n", node->ncvarid); + printf("\tname = %s\n", node->name); + printf("\tnctype = %d\n", node->nctype); + printf("\tndims = %d\n", node->ndims); + printf("\tdimids ="); + for (i = 0; i < node->ndims; i++) + printf(" %d", node->dimids[i]); + printf("\n"); + printf("\tnatts = %d\n", node->natts); + printf("\tmetadata = %d\n", (int)node->metadata); + printf("\tnext = %ld\n", (long)node->next); + printf("\tprev = %ld\n", (long)node->prev); + node = node->next; + } + } + else { + printf("\tList undefined\n"); + } + + printf("End of Variable List\n"); + + return; +} + +void wrap_nc(int status) +{ + if (status != NC_NOERR) { + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); + exit(1); + } + + return; +} + +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len; + struct dim **in_dim_idx; + + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) { + perror("*in_dim_idxp"); + exit(3); + } + in_dim_idx = *in_dim_idxp; + + for (i = 0; i < ndims; i++) { + wrap_nc(nc_inq_dim(ncid, i, name, &len)); + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len); + } + + return; +} +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp) +{ + int i; + char name[NC_MAX_NAME+1], *value; + nc_type xtype; + size_t len; + + for (i = 0; i < ngatts; i++) { + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name)); + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len)); + if (xtype != NC_CHAR) { + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name); + exit(2); + } + if (!(value = (char *)malloc((len+1)*sizeof(char)))) { + perror("get_gatts: value"); + exit(3); + } + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value)); + value[(len+1-1)] = '\0'; + gattlist_add_node(gatt_headp, name, value); + free(value); /* because gattlist_add_node() duplicates it */ + } + + return; +} + +void get_vars(int ncid, int nvars, struct var **var_headp) +{ + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts; + size_t len; + float FillValue; + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag; + nc_type xtype, atype; + + for (i = 0; i < nvars; i++) { + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids, + &natts)); + /* Check for _FillValue */ + FillFlag = 0; + FillValue = 0.; + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len) + == NC_NOERR) { + if (atype == NC_FLOAT && len) { + wrap_nc(nc_get_att_float(ncid, i, + fill_att_name, &FillValue)); + FillFlag = 1; + } + } + varlist_add_node(var_headp, i, name, xtype, ndims, dimids, + natts, FillFlag, FillValue); + } + + return; +} + +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid, + size_t nsamples, int **in2outp, int out_ncid, + struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp) +{ + /* + * Define dimensions on new files and build translation tables between + * dimension IDs. + */ + int j, dimid, ndims, *in2out; + size_t len; + struct dim *dnode; + + ndims = 0; + + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) { + perror("put_dimensions: in2outp"); + exit(3); + } + in2out = *in2outp; + + if (in_dim_head) { + dnode = in_dim_head; + for (j = 0; j == 0 || dnode != in_dim_head; j++) { + if (dnode->dimid == in_unlimdimid) + len = NC_UNLIMITED; + else + len = dnode->len; + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid)); + dimlist_add_node(out_dim_headp, dimid, dnode->name, len); + ++ndims; + in2out[dnode->dimid] = dimid; + if (dnode->dimid == in_unlimdimid) + *out_unlimdimidp = dimid; + /* + * Just after the "time" dimension, add the new + * "nsamples" and "hour" dimensions. + */ + if (!strcmp(dnode->name, time_name)) { + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid)); + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples); + ++ndims; + + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid)); + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY); + ++ndims; + /* Track hour dimid for out files */ + *out_hour_dimidp = dimid; + } + + dnode = dnode->next; + } + } + else { + fprintf(stderr, "WARNING: No dimensions defined!\n"); + } + + return ndims; +} + +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history) +{ + /* + * Write out global attributes matching those of the input file. + * Change history attribute to the string passed in as an argument. + */ + int j, hflag = 0, ngatts; + char *value; + struct text_att *anode; + + ngatts = 0; + + if (in_gatt_head) { + anode = in_gatt_head; + for (j = 0; j == 0 || anode != in_gatt_head; j++) { + if (!strcmp(anode->name, history_name)) { + value = out_history; + ++hflag; + } + else + value = anode->value; + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value)); + ++ngatts; + anode = anode->next; + } + /* If no history attribute on input, add one to the output */ + if (!hflag) { + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history)); + ++ngatts; + } + } + else { + fprintf(stderr, "WARNING: No global attributes defined!\n"); + } + + return ngatts; +} + +int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid) +{ + /* + * Translate between input dimension IDs and those for the output file. + * For normal time-based variables, add a new dimension for hour of + * day. For metadata variables, do not add new dimensions, even if + * it is time-based. Return (possibly new) number of dimensions. + */ + int i, ndims; + + for (i = ndims = 0; i < in_ndims; i++, ndims++) { + out_dimids[ndims] = in2out[in_dimids[i]]; + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name) + && !metadata) { + ndims++; + out_dimids[ndims] = hour_dimid; + } + } + + return ndims; +} + +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid, + int in_varid, int out_ncid, int out_varid) +{ + /* + * Copy the attributes of the input variable to those of the output + * variable. If the variable is not a metadata variable, ensure that + * the cell_method attribute is set either to "time: mean" or + * "time: standard_deviation", depending on the output file type. + */ + + int i, natts; + char name[NC_MAX_NAME+1], cmflag = 0; + char *cell_methods[] = { "time: mean", "time: standard_deviation" }; + + for (i = natts = 0; i < in_natts; i++, natts++) { + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name)); + if (!strcmp(name, cell_method_name) && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + cmflag = 1; + } + else + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid)); + } + /* + * If no cell_method attribute was specified for a non-metadata + * variable on the input file, add the appropriate cell_method anyway + * on the output file. + */ + if (!cmflag && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + natts++; + } + + return natts; +} + +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head, + int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp, + int stat_type) +{ + /* + * Define variables on output file and copy attributes from input file. + * Return number of variables defined. + */ + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts; + struct var *vnode; + + nvars = 0; + + if (in_var_head) { + vnode = in_var_head; + /* + * March through input variables creating (mostly) the same + * variables on the output file. + */ + for (j = 0; j == 0 || vnode != in_var_head; j++) { + ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid); + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid)); + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid); + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + /* + * Just after the "time" variable, add the new "hour" + * variable for hour of day. + */ + if (!strcmp(vnode->name, time_name)) { + ndims = 1; + dimids[0] = hour_dimid; + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid)); + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name)); + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units)); + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0); + ++nvars; + } + + vnode = vnode->next; + } + } + else { + fprintf(stderr, "WARNING: No variables defined!\n"); + } + + return nvars; +} + +void *alloc_var(nc_type nctype, size_t len) +{ + void *val; + + switch (nctype) { + case NC_FLOAT: + if (!(val = (float *)malloc((len) * sizeof(float)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_INT: + if (!(val = (int *)malloc((len) * sizeof(int)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_DOUBLE: + if (!(val = (double *)malloc((len) * sizeof(double)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_CHAR: + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) { + perror("alloc_var: val"); + exit(3); + } + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx) +{ + int i; + size_t len = (size_t)1; + void *val; + + /* Figure out the total size */ + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len; + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_get_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice) +{ + int i; + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + void *val; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + /* + * Figure out the total size for the slice (start at second dimension) + * and build start/count vectors. + */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + len *= (dim_idx[dimids[i]])->len; + } + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void write_var(int ncid, int varid, nc_type nctype, void *val) +{ + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_put_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice) +{ + int i; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + + /* Build start/count vectors */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + } + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void copy_metadata(int in_ncid, struct var *in_var_head, + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head) +{ + int i, j; + float hr[HOURS_PER_DAY]; + struct var *in_vnode, *out_vnode; + void *val; + + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i; + + if (in_var_head) { + in_vnode = in_var_head; + /* + * March through input variables to stuff metadata values into + * the new output files. NOTE: Time-based metadata variables + * should contain only the last (time-slice) value (from all + * input files). + */ + for (j = 0; j == 0 || in_vnode != in_var_head; j++) { + if (in_vnode->metadata) { + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name); + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + /* time is not the first dimension */ +#ifdef DEBUG + printf("Copying metadata variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx); + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val); + } + else { + /* time is the first dimension */ +#ifdef DEBUG + printf("Copying last value of \ +time-based metadata variable %s\n", in_vnode->name); +#endif /* DEBUG */ + 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)); + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0); + } + free(val); + /* + * Just after the "time" variable, write out + * the "hour" variable values. + */ + if (!strcmp(in_vnode->name, time_name)) { + out_vnode = varlist_find_by_name(out_var_head, hour_name); + write_var(out_ncid, out_vnode->ncvarid, + out_vnode->nctype, hr); + } + } + else { +#ifdef DEBUG + printf("Skipping variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + } + in_vnode = in_vnode->next; + } + } + + return; +} + +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples) +{ + char *mean_history_gatt, *stddev_history_gatt, + *mean_prefix="Statistical means from history files:", + *stddev_prefix="Statistical standard deviations from history files:"; + + /* + * Construct strings for history global attributes for the two output + * files. + */ + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:mean_history_gatt"); + exit(2); + } + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:stddev_history_gatt"); + exit(2); + } + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist); + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist); + + /* Open input file */ + wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid)); + /* Inquire about number of dimensions, variables, global attributes + * and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts, + &input_unlimdimid)); + /* Grab dimension IDs and lengths from input file */ + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx); + /* Grab desired global attributes from input file */ + get_gatts(input_ncid, input_ngatts, &input_gatt_head); + /* Grab variable IDs and attributes from input file */ + get_vars(input_ncid, input_nvars, &input_var_head); +#ifdef DEBUG + varlist_print(input_var_head); +#endif /* DEBUG */ + /* Create netCDF files for monthly means and standard deviations */ + /* Create new netCDF files */ + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid)); + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid)); + /* Define dimensions */ + mean_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2mdid, mean_ncid, + &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid); + stddev_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid, + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid); + /* Define variables and their attributes */ + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head, + MEAN_TYPE); + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head, + STDDEV_TYPE); + /* Store global attributes */ + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt); + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid, + stddev_history_gatt); + /* End define mode */ + wrap_nc(nc_enddef(mean_ncid)); + wrap_nc(nc_enddef(stddev_ncid)); + /* Write out metdata variables that do not depend on "time" */ + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid, + mean_var_head); + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid, + stddev_var_head); + + wrap_nc(nc_close(input_ncid)); + + return; +} + +size_t count_samples(int nifnames, char **input_fnames) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len, cnt = 0; + + for (i = 0; i < nifnames; i++) { + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } +#ifdef DEBUG + printf("%ld samples in %s\n", (long)len, input_fnames[i]); +#endif /* DEBUG */ + wrap_nc(nc_close(input_ncid)); + cnt += len; + } + return cnt; +} + +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp) +{ + /* + * Allocate space for arrays of means and standard deviations by + * hour of day. + */ + int i; + size_t **cell_samples; + double **means, **stddevs; + + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*meansp"); + exit(2); + } + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*stddevsp"); + exit(2); + } + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) { + perror("alloc_means_stddevs:*cell_samplesp"); + exit(2); + } + + means = *meansp; + stddevs = *stddevsp; + cell_samples = *cell_samplesp; + + for (i = 0; i < HOURS_PER_DAY; i++) { + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:means[i]"); + exit(2); + } + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:stddevs[i]"); + exit(2); + } + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) { + perror("alloc_means_stddevs:cell_samples[i]"); + exit(2); + } + } + + return; +} + +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue) +{ + int hr, i, j, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + means[hr][pos] = FillValue; + stddevs[hr][pos] = FillValue; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples) +{ + /* + * Free space from arrays of means and standard deviations by + * hour of day. + */ + int i; + + for (i = 0; i < HOURS_PER_DAY; i++) { + free(means[i]); + free(stddevs[i]); + free(cell_samples[i]); + } + + free(means); + free(stddevs); + free(cell_samples); + + return; +} + +void add_to_means_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = (double)val[pos]; + else + means[hr][pos] += (double)val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_means_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = val[pos]; + else + means[hr][pos] += val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + means[hr][pos] /= (double)cell_samples[hr][pos]; + } + } + } + } + + return; +} + +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - (double)val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + stddevs[hr][pos] /= (double)cell_samples[hr][pos]; + stddevs[hr][pos] = sqrt(stddevs[hr][pos]); + } + } + } + } + + return; +} + + +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) +{ + int i, ts; + int varid, *mcsec; + size_t len, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + char name[NC_MAX_NAME+1]; + void *val; + + /* Allocate space for timeslice */ + val = alloc_var(nctype, (d1 * d2)); + + for (i = 0; i < nifnames; i++) { +#ifdef DEBUG + printf("\tOpening %s", input_fnames[i]); + if (ndims > 3) printf(" to retrieve level %d", level); + if (stat_type == MEAN_TYPE) + printf(" for computing means\n"); + else + printf(" for computing stddevs\n"); +#endif /* DEBUG */ + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } + /* Get seconds of day variable */ + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid)); + if (!(mcsec = (int *)malloc(len * sizeof(int)))) { + perror("read_and_compute:mcsec"); + exit(2); + } + wrap_nc(nc_get_var_int(input_ncid, varid, mcsec)); + /* Get variable ID for requested variable */ + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid)); + /* Retrieve time slice of desired variable */ + for (ts = 0; ts < len; ts++) { + start[0] = ts; + count[0] = 1; + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 4) { + start[1] = level; + count[1] = 1; + } + switch(nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, val)); + if (stat_type == MEAN_TYPE) + add_to_means_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_float(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, val)); + if (stat_type == MEAN_TYPE) + add_to_means_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_double(val, mcsec[ts], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + /* Free mcsec vector */ + free(mcsec); + /* Close input file */ + wrap_nc(nc_close(input_ncid)); + } + /* Divide sums by number of samples */ + if (stat_type == MEAN_TYPE) + divide_means(d1, d2, means, cell_samples); + else + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples); + + /* Free working variable array */ + free(val); + + return; +} + +float *double_to_float(size_t len, double *dvar) +{ + int i; + float *fvar; + + if (!(fvar = (float *)malloc(len * sizeof(float)))) { + perror("double_to_float:fvar"); + exit(2); + } + + for (i = 0; i < len; i++) + fvar[i] = (float)dvar[i]; + + return fvar; +} + +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2, + int level, int ndims, double **var_hours) +{ + /* Output dimensions should be one larger than input dimensions */ + int i, hr; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + float *fvar = NULL; + + if (nctype == NC_FLOAT) { + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) { + perror("write_var_hours:fvar"); + exit(2); + } + } + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + start[0] = 0; + count[0] = 1; /* time */ + start[1] = hr; + count[1] = 1; /* hour */ + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 5) { + start[2] = level; + count[2] = 1; + } + switch (nctype) { + case NC_FLOAT: + for (i = 0; i < (d1 * d2); i++) + fvar[i] = (float)var_hours[hr][i]; + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr])); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + if (nctype == NC_FLOAT) + free(fvar); + + return; +} + +void compute_stats(int nifnames, char **input_fnames, size_t nsamples) +{ + int j, k, nlevels; + size_t d1, d2, **cell_samples; + double **means, **stddevs; + struct var *in_vnode, *out_vnode; + + if (input_var_head) { + in_vnode = input_var_head; + /* March through non-metadata variables to compute stats */ + for (j = 0; j == 0 || in_vnode != input_var_head; j++) { + if (!in_vnode->metadata) { + /* Make sure time is really the first dimension */ + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name); + exit(3); + } + /* Figure out input dimensions */ + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) { + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims); + exit(3); + } + /* Assume 2D output; find dimensions */ + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len; + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len; + if (in_vnode->ndims == 3) + nlevels = 1; + else + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len; + /* Allocate working space for means and stddevs */ + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples); + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue); + printf("Computing statistics for %s\n", + in_vnode->name); + /* Compute means and stddevs, then write them */ + for (k = 0; k < nlevels; k++) { + /* Read and compute means */ + 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); + /* Find corresponding output variable on the mean output file */ + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name); + /* Write out the means for this variable */ + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means); + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */ + reset_cell_samples(d1, d2, cell_samples); + /* Read and compute stddevs using means */ + 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); + /* Find corresponding output variable on the stddev output file */ + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name); + /* Write out stddevs for this variable */ + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs); + } + + /* Free working space for means and stddevs */ + free_means_stddevs(means, stddevs, cell_samples); + } + in_vnode = in_vnode->next; + } + } + return; +} + +void usage(char *program) +{ + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program); + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n"); + exit(4); +} + +int main(int argc, char **argv) +{ + int i, ic, nifnames; + size_t len, pos, nsamples; + char *mfname, *sfname, **ifnames, *flist; + + ifnames = NULL; + mfname = sfname = flist = NULL; + nifnames = 0; + +#ifdef DEBUG + printf("Number of metadata variables (nmvars) = %d\n", nmvars); + fflush(stdout); +#endif /* DEBUG */ + + /* check command-line flags and switches */ + while ((ic = getopt(argc, argv, "m:s:")) != -1) { + switch(ic) { + case 'm': + if (!(mfname = strdup(optarg))) { + perror("mfname"); + exit(2); + } + break; + case 's': + if (!(sfname = strdup(optarg))) { + perror("sfname"); + exit(2); + } + break; + } + } + + if (!mfname) { + fprintf(stderr, "Output file name for writing means required.\n"); + usage(argv[0]); + } + if (!sfname) { + fprintf(stderr, "Output file name for writing standard deviations required.\n"); + usage(argv[0]); + } + + if (optind < argc) { + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) { + perror("ifnames"); + exit(2); + } + for (i = optind; i < argc; i++) { + if (!(ifnames[nifnames++] = strdup(argv[i]))) { + perror("ifnames[nifnames++]"); + exit(2); + } + } + } + else { + fprintf(stderr, "At least one input file name is required.\n"); + usage(argv[0]); + } + + + /* + * Build list of input files to be included in the global history + * attribute on the two outputs files. + */ + for (i = len = 0; i < nifnames; i++) + len += strlen(ifnames[i]); + len += nifnames + 1; /* space in front of every name + null terminator */ + if (!(flist = (char *)malloc(len * sizeof(char)))) { + perror("flist"); + exit(2); + } + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++) + sprintf(&flist[pos], " %s", ifnames[i]); +#ifdef DEBUG + printf("flist=%s\n", flist); + fflush(stdout); +#endif /* DEBUG */ + + /* Open every file to count up the number of time samples. */ + nsamples = count_samples(nifnames, ifnames); +#ifdef DEBUG + printf("Number of samples across all files = %ld\n", (long)nsamples); + fflush(stdout); +#endif /* DEBUG */ + + /* + * Open the *last* input file and the two output files (for mean and + * standard deviation). Define dimensions, variables, and attributes + * for the two output files based on the *last* input files. The + * last file is used because some metadata variables will contain + * only values from the last time slice from the period over which + * the statistics are computed. + */ + open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples); + + compute_stats(nifnames, ifnames, nsamples); + + /* Close the two output files */ + wrap_nc(nc_close(mean_ncid)); + wrap_nc(nc_close(stddev_ncid)); + + return 0; +} diff -r 000000000000 -r 3c02cce30be8 h1_summary2.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/h1_summary2.c Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,1614 @@ +#include +#include +#include +#include +#include +#include "netcdf.h" + +/* + * Loop through one month of hourly history tapes from the Community Land Model + * and generate monthly statistics (means and standard deviations) of fields + * by hour of day. + */ + +#define MEAN_TYPE 0 +#define STDDEV_TYPE 1 + +#define HOURS_PER_DAY 24 +#define MIN_PER_HOUR 60 +#define SEC_PER_MIN 60 +#define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN) + +static char *metadata_vars[] = { + "lon", + "lat", + "lonrof", + "latrof", + "time", + "hour", /* new metadata variable to be added to output files */ + "levsoi", + "levlak", + "edgen", + "edgee", + "edges", + "edgew", + "longxy", + "latixy", + "area", + "areaupsc", + "topo", + "topodnsc", + "landfrac", + "landmask", + "pftmask", + "indxupsc", + "mcdate", + "mcsec", + "mdcur", + "mscur", + "nstep", + "time_bounds", + "date_written", + "time_written" +}; + +struct text_att { + char *name; + char *value; + struct text_att *next; + struct text_att *prev; +}; + +struct dim { + int dimid; + char *name; + size_t len; + struct dim *next; + struct dim *prev; +}; + +struct var { + int ncvarid; + char *name; + nc_type nctype; + int ndims; + int *dimids; + int natts; + char metadata; + char FillFlag; + float FillValue; + struct var *next; + struct var *prev; +}; + +static char *time_name = "time"; +static char *mcsec_name = "mcsec"; +static char *history_name = "history"; +static char *nsamples_name = "nsamples"; +static char *hour_name = "hour", *hour_long_name = "hour of day", + *hour_units = "hour"; +static char *cell_method_name = "cell_method"; +/* input stuff */ +static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars); +static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid; +static struct text_att *input_gatt_head = NULL; +static struct dim *input_dim_head = NULL, **input_dim_idx; +static struct var *input_var_head = NULL; +/* translation stuff */ +static int *idid2mdid, *idid2sdid; /* dimension IDs */ +/* output stuff */ +static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid; +static int mean_hour_dimid; /* special notes */ +static struct dim *mean_dim_head = NULL; +static struct var *mean_var_head = NULL; +static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid; +static int stddev_hour_dimid; /* special notes */ +static struct dim *stddev_dim_head = NULL; +static struct var *stddev_var_head = NULL; + +char is_metadata(char *name) +{ + int i; + + for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++); + + if (i < nmvars) + return 1; + else + return 0; +} + +#if 0 +struct dim *dimlist_find_by_name(struct dim *head, char *name) +{ + int i; + struct dim *p = NULL; + + if (head) { + node = head; + for (i = 0 ; i == 0 || node != head; i++) { + if (!strcmp(node->name, name)) + p = node; + node = node->next; + } + return p; + } + else + return NULL; +} +#endif + +struct var *varlist_find_by_name(struct var *head, char *name) +{ + int i; + struct var *node; + + if (head) { + node = head; + for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++) + node = node->next; + if (i && node == head) + return NULL; + else + return node; + } + else + return NULL; +} + +void gattlist_add_node(struct text_att **headp, char *name, char *value) +{ + struct text_att *head, *node; + + head = *headp; + + if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) { + perror("gattlist_add_node"); + exit(2); + } + + if (!(node->name = strdup(name))) { + perror("gattlist_add_node"); + exit(2); + } + if (!(node->value = strdup(value))) { + perror("gattlist_add_node"); + exit(2); + } + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len) +{ + struct dim *head, *node; + + head = *headp; + + if (!(node = (struct dim *)malloc(sizeof(struct dim)))) { + perror("dimlist_add_node"); + exit(2); + } + + node->dimid = dimid; + if (!(node->name = strdup(name))) { + perror("dimlist_add_node"); + exit(2); + } + node->len = len; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return node; +} + +void varlist_add_node(struct var **headp, int ncvarid, char *name, + nc_type nctype, int ndims, int *dimids, int natts, char FillFlag, + float FillValue) +{ + int i; + struct var *head, *node; + + head = *headp; + + if (!(node = (struct var *)malloc(sizeof(struct var)))) { + perror("varlist_add_node"); + exit(3); + } + + node->ncvarid = ncvarid; + if (!(node->name = strdup(name))) { + perror("varlist_add_node"); + exit(3); + } + node->nctype = nctype; + node->ndims = ndims; + if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) { + perror("varlist_add_node"); + exit(3); + } + for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i]; + node->natts = natts; + node->metadata = is_metadata(name); + node->FillFlag = FillFlag; + node->FillValue = FillValue; + + if (head == NULL) { + /* first one */ + *headp = node; + node->next = node; + node->prev = node; + } + else { + node->next = head; + node->prev = head->prev; + /* set this after setting node->prev from here */ + head->prev = node; + /* set this after having set node->prev */ + node->prev->next = node; + } + + return; +} + +void varlist_print(struct var *headp) +{ + int i, j; + struct var *node; + + printf("Beginning of Variable List\n"); + + if (headp) { + node = headp; + for (j = 0; j == 0 || node != headp; j++) { + printf("Variable %d (ptr=%ld):\n", j, (long)node); + printf("\tncvarid = %d\n", node->ncvarid); + printf("\tname = %s\n", node->name); + printf("\tnctype = %d\n", node->nctype); + printf("\tndims = %d\n", node->ndims); + printf("\tdimids ="); + for (i = 0; i < node->ndims; i++) + printf(" %d", node->dimids[i]); + printf("\n"); + printf("\tnatts = %d\n", node->natts); + printf("\tmetadata = %d\n", (int)node->metadata); + printf("\tnext = %ld\n", (long)node->next); + printf("\tprev = %ld\n", (long)node->prev); + node = node->next; + } + } + else { + printf("\tList undefined\n"); + } + + printf("End of Variable List\n"); + + return; +} + +void wrap_nc(int status) +{ + if (status != NC_NOERR) { + fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); + exit(1); + } + + return; +} + +void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len; + struct dim **in_dim_idx; + + if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) { + perror("*in_dim_idxp"); + exit(3); + } + in_dim_idx = *in_dim_idxp; + + for (i = 0; i < ndims; i++) { + wrap_nc(nc_inq_dim(ncid, i, name, &len)); + in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len); + } + + return; +} +void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp) +{ + int i; + char name[NC_MAX_NAME+1], *value; + nc_type xtype; + size_t len; + + for (i = 0; i < ngatts; i++) { + wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name)); + wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len)); + if (xtype != NC_CHAR) { + fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name); + exit(2); + } + if (!(value = (char *)malloc((len+1)*sizeof(char)))) { + perror("get_gatts: value"); + exit(3); + } + wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value)); + value[(len+1-1)] = '\0'; + gattlist_add_node(gatt_headp, name, value); + free(value); /* because gattlist_add_node() duplicates it */ + } + + return; +} + +void get_vars(int ncid, int nvars, struct var **var_headp) +{ + int i, ndims, dimids[NC_MAX_VAR_DIMS], natts; + size_t len; + float FillValue; + char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag; + nc_type xtype, atype; + + for (i = 0; i < nvars; i++) { + wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids, + &natts)); + /* Check for _FillValue */ + FillFlag = 0; + FillValue = 0.; + if (nc_inq_att(ncid, i, fill_att_name, &atype, &len) + == NC_NOERR) { + if (atype == NC_FLOAT && len) { + wrap_nc(nc_get_att_float(ncid, i, + fill_att_name, &FillValue)); + FillFlag = 1; + } + } + varlist_add_node(var_headp, i, name, xtype, ndims, dimids, + natts, FillFlag, FillValue); + } + + return; +} + +int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid, + size_t nsamples, int **in2outp, int out_ncid, + struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp) +{ + /* + * Define dimensions on new files and build translation tables between + * dimension IDs. + */ + int j, dimid, ndims, *in2out; + size_t len; + struct dim *dnode; + + ndims = 0; + + if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) { + perror("put_dimensions: in2outp"); + exit(3); + } + in2out = *in2outp; + + if (in_dim_head) { + dnode = in_dim_head; + for (j = 0; j == 0 || dnode != in_dim_head; j++) { + if (dnode->dimid == in_unlimdimid) + len = NC_UNLIMITED; + else + len = dnode->len; + wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid)); + dimlist_add_node(out_dim_headp, dimid, dnode->name, len); + ++ndims; + in2out[dnode->dimid] = dimid; + if (dnode->dimid == in_unlimdimid) + *out_unlimdimidp = dimid; + /* + * Just after the "time" dimension, add the new + * "nsamples" and "hour" dimensions. + */ + if (!strcmp(dnode->name, time_name)) { + wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid)); + dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples); + ++ndims; + + wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid)); + dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY); + ++ndims; + /* Track hour dimid for out files */ + *out_hour_dimidp = dimid; + } + + dnode = dnode->next; + } + } + else { + fprintf(stderr, "WARNING: No dimensions defined!\n"); + } + + return ndims; +} + +int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history) +{ + /* + * Write out global attributes matching those of the input file. + * Change history attribute to the string passed in as an argument. + */ + int j, hflag = 0, ngatts; + char *value; + struct text_att *anode; + + ngatts = 0; + + if (in_gatt_head) { + anode = in_gatt_head; + for (j = 0; j == 0 || anode != in_gatt_head; j++) { + if (!strcmp(anode->name, history_name)) { + value = out_history; + ++hflag; + } + else + value = anode->value; + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value)); + ++ngatts; + anode = anode->next; + } + /* If no history attribute on input, add one to the output */ + if (!hflag) { + wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history)); + ++ngatts; + } + } + else { + fprintf(stderr, "WARNING: No global attributes defined!\n"); + } + + return ngatts; +} + +int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid) +{ + /* + * Translate between input dimension IDs and those for the output file. + * For normal time-based variables, add a new dimension for hour of + * day. For metadata variables, do not add new dimensions, even if + * it is time-based. Return (possibly new) number of dimensions. + */ + int i, ndims; + + for (i = ndims = 0; i < in_ndims; i++, ndims++) { + out_dimids[ndims] = in2out[in_dimids[i]]; + if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name) + && !metadata) { + ndims++; + out_dimids[ndims] = hour_dimid; + } + } + + return ndims; +} + +int copy_att(char metadata, int stat_type, int in_natts, int in_ncid, + int in_varid, int out_ncid, int out_varid) +{ + /* + * Copy the attributes of the input variable to those of the output + * variable. If the variable is not a metadata variable, ensure that + * the cell_method attribute is set either to "time: mean" or + * "time: standard_deviation", depending on the output file type. + */ + + int i, natts; + char name[NC_MAX_NAME+1], cmflag = 0; + char *cell_methods[] = { "time: mean", "time: standard_deviation" }; + + for (i = natts = 0; i < in_natts; i++, natts++) { + wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name)); + if (!strcmp(name, cell_method_name) && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + cmflag = 1; + } + else + wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid)); + } + /* + * If no cell_method attribute was specified for a non-metadata + * variable on the input file, add the appropriate cell_method anyway + * on the output file. + */ + if (!cmflag && !metadata) { + wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); + natts++; + } + + return natts; +} + +int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head, + int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp, + int stat_type) +{ + /* + * Define variables on output file and copy attributes from input file. + * Return number of variables defined. + */ + int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts; + struct var *vnode; + + nvars = 0; + + if (in_var_head) { + vnode = in_var_head; + /* + * March through input variables creating (mostly) the same + * variables on the output file. + */ + for (j = 0; j == 0 || vnode != in_var_head; j++) { + ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid); + wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid)); + natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid); + varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); + ++nvars; + /* + * Just after the "time" variable, add the new "hour" + * variable for hour of day. + */ + if (!strcmp(vnode->name, time_name)) { + ndims = 1; + dimids[0] = hour_dimid; + wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid)); + wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name)); + wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units)); + varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0); + ++nvars; + } + + vnode = vnode->next; + } + } + else { + fprintf(stderr, "WARNING: No variables defined!\n"); + } + + return nvars; +} + +void *alloc_var(nc_type nctype, size_t len) +{ + void *val; + + switch (nctype) { + case NC_FLOAT: + if (!(val = (float *)malloc((len) * sizeof(float)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_INT: + if (!(val = (int *)malloc((len) * sizeof(int)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_DOUBLE: + if (!(val = (double *)malloc((len) * sizeof(double)))) { + perror("alloc_var: val"); + exit(3); + } + break; + case NC_CHAR: + if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) { + perror("alloc_var: val"); + exit(3); + } + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx) +{ + int i; + size_t len = (size_t)1; + void *val; + + /* Figure out the total size */ + for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len; + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_get_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice) +{ + int i; + size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + void *val; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + /* + * Figure out the total size for the slice (start at second dimension) + * and build start/count vectors. + */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + len *= (dim_idx[dimids[i]])->len; + } + + val = alloc_var(nctype, len); + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_get_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_get_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return val; +} + +void write_var(int ncid, int varid, nc_type nctype, void *val) +{ + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_var_float(ncid, varid, val)); + break; + case NC_INT: + wrap_nc(nc_put_var_int(ncid, varid, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_var_double(ncid, varid, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_var_text(ncid, varid, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice) +{ + int i; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + + /* Make sure time is really the first dimension */ + if (strcmp((dim_idx[dimids[0]])->name, time_name)) { + fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name); + exit(3); + } + + /* Build start/count vectors */ + start[0] = tslice; + count[0] = 1; + for (i = 1; i < ndims; i++) { + start[i] = 0; + count[i] = (dim_idx[dimids[i]])->len; + } + + switch (nctype) { + case NC_FLOAT: + wrap_nc(nc_put_vara_float(ncid, varid, start, count, val)); + break; + case NC_INT: + wrap_nc(nc_put_vara_int(ncid, varid, start, count, val)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, val)); + break; + case NC_CHAR: + wrap_nc(nc_put_vara_text(ncid, varid, start, count, val)); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + return; +} + +void copy_metadata(int in_ncid, struct var *in_var_head, + struct dim **in_dim_idx, int out_ncid, struct var *out_var_head) +{ + int i, j; + float hr[HOURS_PER_DAY]; + struct var *in_vnode, *out_vnode; + void *val; + + for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i; + + if (in_var_head) { + in_vnode = in_var_head; + /* + * March through input variables to stuff metadata values into + * the new output files. NOTE: Time-based metadata variables + * should contain only the last (time-slice) value (from all + * input files). + */ + for (j = 0; j == 0 || in_vnode != in_var_head; j++) { + if (in_vnode->metadata) { + out_vnode = varlist_find_by_name(out_var_head, in_vnode->name); + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + /* time is not the first dimension */ +#ifdef DEBUG + printf("Copying metadata variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx); + write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val); + } + else { + /* time is the first dimension */ +#ifdef DEBUG + printf("Copying last value of \ +time-based metadata variable %s\n", in_vnode->name); +#endif /* DEBUG */ + 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)); + write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0); + } + free(val); + /* + * Just after the "time" variable, write out + * the "hour" variable values. + */ + if (!strcmp(in_vnode->name, time_name)) { + out_vnode = varlist_find_by_name(out_var_head, hour_name); + write_var(out_ncid, out_vnode->ncvarid, + out_vnode->nctype, hr); + } + } + else { +#ifdef DEBUG + printf("Skipping variable %s\n", + in_vnode->name); +#endif /* DEBUG */ + } + in_vnode = in_vnode->next; + } + } + + return; +} + +void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples) +{ + char *mean_history_gatt, *stddev_history_gatt, + *mean_prefix="Statistical means from history files:", + *stddev_prefix="Statistical standard deviations from history files:"; + + /* + * Construct strings for history global attributes for the two output + * files. + */ + if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:mean_history_gatt"); + exit(2); + } + if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) { + perror("open_inout:stddev_history_gatt"); + exit(2); + } + sprintf(mean_history_gatt, "%s%s", mean_prefix, flist); + sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist); + + /* Open input file */ + wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid)); + /* Inquire about number of dimensions, variables, global attributes + * and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts, + &input_unlimdimid)); + /* Grab dimension IDs and lengths from input file */ + get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx); + /* Grab desired global attributes from input file */ + get_gatts(input_ncid, input_ngatts, &input_gatt_head); + /* Grab variable IDs and attributes from input file */ + get_vars(input_ncid, input_nvars, &input_var_head); +#ifdef DEBUG + varlist_print(input_var_head); +#endif /* DEBUG */ + /* Create netCDF files for monthly means and standard deviations */ + /* Create new netCDF files */ + wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid)); + wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid)); + /* Define dimensions */ + mean_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2mdid, mean_ncid, + &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid); + stddev_ndims = put_dimensions(input_dim_head, input_ndims, + input_unlimdimid, nsamples, &idid2sdid, stddev_ncid, + &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid); + /* Define variables and their attributes */ + mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head, + MEAN_TYPE); + stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, + idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head, + STDDEV_TYPE); + /* Store global attributes */ + mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt); + stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid, + stddev_history_gatt); + /* End define mode */ + wrap_nc(nc_enddef(mean_ncid)); + wrap_nc(nc_enddef(stddev_ncid)); + /* Write out metdata variables that do not depend on "time" */ + copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid, + mean_var_head); + copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid, + stddev_var_head); + + wrap_nc(nc_close(input_ncid)); + + return; +} + +size_t count_samples(int nifnames, char **input_fnames) +{ + int i; + char name[NC_MAX_NAME+1]; + size_t len, cnt = 0; + + for (i = 0; i < nifnames; i++) { + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } +#ifdef DEBUG + printf("%ld samples in %s\n", (long)len, input_fnames[i]); +#endif /* DEBUG */ + wrap_nc(nc_close(input_ncid)); + cnt += len; + } + return cnt; +} + +void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp) +{ + /* + * Allocate space for arrays of means and standard deviations by + * hour of day. + */ + int i; + size_t **cell_samples; + double **means, **stddevs; + + if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*meansp"); + exit(2); + } + if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { + perror("alloc_means_stddevs:*stddevsp"); + exit(2); + } + if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) { + perror("alloc_means_stddevs:*cell_samplesp"); + exit(2); + } + + means = *meansp; + stddevs = *stddevsp; + cell_samples = *cell_samplesp; + + for (i = 0; i < HOURS_PER_DAY; i++) { + if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:means[i]"); + exit(2); + } + if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { + perror("alloc_means_stddevs:stddevs[i]"); + exit(2); + } + if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) { + perror("alloc_means_stddevs:cell_samples[i]"); + exit(2); + } + } + + return; +} + +void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue) +{ + int hr, i, j, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + means[hr][pos] = FillValue; + stddevs[hr][pos] = FillValue; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + cell_samples[hr][pos] = 0; + } + } + } + + return; +} + +void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples) +{ + /* + * Free space from arrays of means and standard deviations by + * hour of day. + */ + int i; + + for (i = 0; i < HOURS_PER_DAY; i++) { + free(means[i]); + free(stddevs[i]); + free(cell_samples[i]); + } + + free(means); + free(stddevs); + free(cell_samples); + + return; +} + +void add_to_means_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = (double)val[pos]; + else + means[hr][pos] += (double)val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_means_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue)) { + if (cell_samples[hr][pos] == 0) + means[hr][pos] = val[pos]; + else + means[hr][pos] += val[pos]; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + means[hr][pos] /= (double)cell_samples[hr][pos]; + } + } + } + } + + return; +} + +void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - (double)val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + +void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2, + char FillFlag, float FillValue, double **means, double **stddevs, + size_t **cell_samples) +{ + int i, j, hr, pos; + double delta; + + hr = (int)floor((double)sec/(double)SEC_PER_HOUR); + + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (!FillFlag || (FillFlag && val[pos] != FillValue + && means[hr][pos] != FillValue)) { + delta = means[hr][pos] - val[pos]; + if (cell_samples[hr][pos] == 0) + stddevs[hr][pos] = delta * delta; + else + stddevs[hr][pos] += delta * delta; + ++cell_samples[hr][pos]; + } + } + } + + return; +} + + +void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples) +{ + int i, j, hr, pos; + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + for (i = 0; i < d1; i++) { +#pragma _CRI concurrent + for (j = 0; j < d2; j++) { + pos = i*d2+j; + if (cell_samples[hr][pos] != 0) { + stddevs[hr][pos] /= (double)cell_samples[hr][pos]; + stddevs[hr][pos] = sqrt(stddevs[hr][pos]); + } + } + } + } + + return; +} + +void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp) +{ + int i; + void **vals; + + /* Allocate space for all timeslices */ + if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) { + perror("alloc_all_samples:*valsp"); + exit(2); + } + vals = *valsp; + for (i = 0; i < nsamples; i++) { + vals[i] = alloc_var(nctype, (d1 * d2)); + } + if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) { + perror("alloc_all_samples:*mcsecp"); + exit(2); + } + + return; +} + +void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec) +{ + int i, j; + float **fvals = NULL; + double **dvals = NULL; + + switch(nctype) { + case NC_FLOAT: + fvals = (float **)vals; + break; + case NC_DOUBLE: + dvals = (double **)vals; + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + + for (i = 0; i < nsamples; i++) { + for (j = 0; j < (d1 * d2); j++) { + switch(nctype) { + case NC_FLOAT: + fvals[i][j] = (float)0.0; + break; + case NC_DOUBLE: + dvals[i][j] = (double)0.0; + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + mcsec[i] = 0; + } + + return; +} + +void free_all_samples(size_t nsamples, void **vals, int *mcsec) +{ + int i; + + for (i = 0; i < nsamples; i++) + if (vals[i]) free(vals[i]); + + free(vals); + free(mcsec); + + return; +} + +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) +{ + int i, ts; + int varid; + size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + char name[NC_MAX_NAME+1]; + + for (i = 0; i < nifnames; i++) { +#ifdef DEBUG + printf("\tOpening %s", input_fnames[i]); + if (ndims > 3) printf(" to retrieve level %d\n", level); +#endif /* DEBUG */ + /* Open input file */ + wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); + /* + * Inquire about number of dimensions, variables, global + * attributes and the ID of the unlimited dimension + */ + wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, + &input_ngatts, &input_unlimdimid)); + wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); + if (strcmp(name, time_name)) { + fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); + exit(4); + } + /* Make sure we don't run off the end of allocated space */ + if ((nslice+len) > nsamples) { + fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples); + exit(4); + } + /* Get seconds of day variable */ + wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid)); + wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice])); + /* Get variable ID for requested variable */ + wrap_nc(nc_inq_varid(input_ncid, var_name, &varid)); + /* Retrieve time slice of desired variable */ + for (ts = 0; ts < len; ts++) { + start[0] = ts; + count[0] = 1; + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 4) { + start[1] = level; + count[1] = 1; + } + switch(nctype) { + case NC_FLOAT: + wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice])); + break; + case NC_DOUBLE: + wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice])); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + nslice++; + } + + /* Close input file */ + wrap_nc(nc_close(input_ncid)); + } + + return; +} + +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) +{ + int i; + + for (i = 0; i < nsamples; i++) { + switch(nctype) { + case NC_FLOAT: + if (stat_type == MEAN_TYPE) + add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + case NC_DOUBLE: + if (stat_type == MEAN_TYPE) + add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples); + else + add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + /* Divide sums by number of samples */ + if (stat_type == MEAN_TYPE) + divide_means(d1, d2, means, cell_samples); + else + divide_sqrt_stddevs(d1, d2, stddevs, cell_samples); + + return; +} + +float *double_to_float(size_t len, double *dvar) +{ + int i; + float *fvar; + + if (!(fvar = (float *)malloc(len * sizeof(float)))) { + perror("double_to_float:fvar"); + exit(2); + } + + for (i = 0; i < len; i++) + fvar[i] = (float)dvar[i]; + + return fvar; +} + +void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2, + int level, int ndims, double **var_hours) +{ + /* Output dimensions should be one larger than input dimensions */ + int i, hr; + size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; + float *fvar = NULL; + + if (nctype == NC_FLOAT) { + if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) { + perror("write_var_hours:fvar"); + exit(2); + } + } + + for (hr = 0; hr < HOURS_PER_DAY; hr++) { + start[0] = 0; + count[0] = 1; /* time */ + start[1] = hr; + count[1] = 1; /* hour */ + start[(ndims-2)] = 0; + count[(ndims-2)] = d1; + start[(ndims-1)] = 0; + count[(ndims-1)] = d2; + if (ndims == 5) { + start[2] = level; + count[2] = 1; + } + switch (nctype) { + case NC_FLOAT: + for (i = 0; i < (d1 * d2); i++) + fvar[i] = (float)var_hours[hr][i]; + wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar)); + break; + case NC_DOUBLE: + wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr])); + break; + default: + fprintf(stderr, "netCDF external data type %d not supported\n", nctype); + exit(3); + } + } + + if (nctype == NC_FLOAT) + free(fvar); + + return; +} + +void compute_stats(int nifnames, char **input_fnames, size_t nsamples) +{ + int j, k, nlevels, *mcsec; + size_t d1, d2, **cell_samples; + double **means, **stddevs; + struct var *in_vnode, *out_vnode; + void **vals; + + if (input_var_head) { + in_vnode = input_var_head; + /* March through non-metadata variables to compute stats */ + for (j = 0; j == 0 || in_vnode != input_var_head; j++) { + if (!in_vnode->metadata) { + /* Make sure time is really the first dimension */ + if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { + fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name); + exit(3); + } + /* Figure out input dimensions */ + if (in_vnode->ndims < 3 || in_vnode->ndims > 4) { + fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims); + exit(3); + } + /* Assume 2D output; find dimensions */ + d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len; + d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len; + if (in_vnode->ndims == 3) + nlevels = 1; + else + nlevels = (input_dim_idx[in_vnode->dimids[1]])->len; + /* Allocate working space for means and stddevs */ + alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples); + init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue); + /* Allocate and initize space for entire field across all files */ + alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec); + init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec); + printf("Computing statistics for %s\n", + in_vnode->name); + /* Compute means and stddevs, then write them */ + for (k = 0; k < nlevels; k++) { + /* Read the entire fields from all files */ + read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec); + /* Compute means */ + compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples); + /* Find corresponding output variable on the mean output file */ + out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name); + /* Write out the means for this variable */ + write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means); + /* Zero out cell_samples so they can be used as a flag again for computing stddevs */ + reset_cell_samples(d1, d2, cell_samples); + /* Compute stddevs using means */ + compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples); + /* Find corresponding output variable on the stddev output file */ + out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name); + /* Write out stddevs for this variable */ + write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs); + } + + /* Free all samples */ + free_all_samples(nsamples, vals, mcsec); + /* Free working space for means and stddevs */ + free_means_stddevs(means, stddevs, cell_samples); + } + in_vnode = in_vnode->next; + } + } + return; +} + +void usage(char *program) +{ + fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program); + fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n"); + exit(4); +} + +int main(int argc, char **argv) +{ + int i, ic, nifnames; + size_t len, pos, nsamples; + char *mfname, *sfname, **ifnames, *flist; + + ifnames = NULL; + mfname = sfname = flist = NULL; + nifnames = 0; + +#ifdef DEBUG + printf("Number of metadata variables (nmvars) = %d\n", nmvars); + fflush(stdout); +#endif /* DEBUG */ + + /* check command-line flags and switches */ + while ((ic = getopt(argc, argv, "m:s:")) != -1) { + switch(ic) { + case 'm': + if (!(mfname = strdup(optarg))) { + perror("mfname"); + exit(2); + } + break; + case 's': + if (!(sfname = strdup(optarg))) { + perror("sfname"); + exit(2); + } + break; + } + } + + if (!mfname) { + fprintf(stderr, "Output file name for writing means required.\n"); + usage(argv[0]); + } + if (!sfname) { + fprintf(stderr, "Output file name for writing standard deviations required.\n"); + usage(argv[0]); + } + + if (optind < argc) { + if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) { + perror("ifnames"); + exit(2); + } + for (i = optind; i < argc; i++) { + if (!(ifnames[nifnames++] = strdup(argv[i]))) { + perror("ifnames[nifnames++]"); + exit(2); + } + } + } + else { + fprintf(stderr, "At least one input file name is required.\n"); + usage(argv[0]); + } + + + /* + * Build list of input files to be included in the global history + * attribute on the two outputs files. + */ + for (i = len = 0; i < nifnames; i++) + len += strlen(ifnames[i]); + len += nifnames + 1; /* space in front of every name + null terminator */ + if (!(flist = (char *)malloc(len * sizeof(char)))) { + perror("flist"); + exit(2); + } + for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++) + sprintf(&flist[pos], " %s", ifnames[i]); +#ifdef DEBUG + printf("flist=%s\n", flist); + fflush(stdout); +#endif /* DEBUG */ + + /* Open every file to count up the number of time samples. */ + nsamples = count_samples(nifnames, ifnames); +#ifdef DEBUG + printf("Number of samples across all files = %ld\n", (long)nsamples); + fflush(stdout); +#endif /* DEBUG */ + + /* + * Open the *last* input file and the two output files (for mean and + * standard deviation). Define dimensions, variables, and attributes + * for the two output files based on the *last* input files. The + * last file is used because some metadata variables will contain + * only values from the last time slice from the period over which + * the statistics are computed. + */ + open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples); + + compute_stats(nifnames, ifnames, nsamples); + + /* Close the two output files */ + wrap_nc(nc_close(mean_ncid)); + wrap_nc(nc_close(stddev_ncid)); + + return 0; +} diff -r 000000000000 -r 3c02cce30be8 setup_phoenix.bash --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/setup_phoenix.bash Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,17 @@ +#!/bin/bash +source /opt/modules/modules/init/bash +# for robin +#module unload PrgEnv +##module load pgi/6.1 +##module load netcdf_pgi/3.6.1 +#module load netcdf_gcc/3.6.1 +# for phoenix +module purge +module load open +module load PrgEnv.5407 +module unload mpt +module load mpt.2.4.0.6 +module load pbs +module load netcdf/3.5.1_r4 +# +module list diff -r 000000000000 -r 3c02cce30be8 setup_robin1.bash --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/setup_robin1.bash Wed Sep 26 17:16:40 2007 -0400 @@ -0,0 +1,11 @@ +#!/bin/bash +source /opt/modules/modules/init/bash +# for robin1 +module unload PrgEnv +#module load pgi/6.1 +#module load netcdf_pgi/3.6.1 +module load netcdf_gcc/3.6.1 +# for phoenix +#module load netcdf/3.5.1_r4 +# +module list