add_total_fields.c
changeset 0 3c02cce30be8
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/add_total_fields.c	Wed Sep 26 17:16:40 2007 -0400
     1.3 @@ -0,0 +1,548 @@
     1.4 +#include <stdio.h>
     1.5 +#include <stdlib.h>
     1.6 +#include <string.h>
     1.7 +#include <math.h>
     1.8 +#include "netcdf.h"
     1.9 +
    1.10 +/*
    1.11 + * Loop through all history tapes from the Community Land Model adding new
    1.12 + * fields that are totals of the multi-level fields included in total_vars[]
    1.13 + * below.
    1.14 + */
    1.15 +
    1.16 +static char *total_vars[] = {
    1.17 +	"SOILICE",
    1.18 +	"SOILLIQ"
    1.19 +};
    1.20 +static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars);
    1.21 +static char *combo_vars[] = {
    1.22 +	"LATENT",		/* FCTR + FCEV + FGEV */
    1.23 +	"ET",			/* QVEGE + QVEGT + QSOIL */
    1.24 +	"ALL_SKY_ALBEDO",	/* FSR / FSDS */
    1.25 +	"BLACK_SKY_ALBEDO",	/* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */
    1.26 +	"NETRAD",		/* FSA - FIRA */
    1.27 +};
    1.28 +static char *combo_long_names[] = {
    1.29 +	"latent heat flux",
    1.30 +	"evapotranspiration",
    1.31 +	"surface all-sky albedo",
    1.32 +	"surface black-sky albedo",
    1.33 +	"net radiation"
    1.34 +};
    1.35 +static char *combo_units[] = {
    1.36 +	"watt/m^2",
    1.37 +	"m/s",
    1.38 +	"unitless",
    1.39 +	"unitless",
    1.40 +	"watt/m^2"
    1.41 +};
    1.42 +static float combo_FillValues[] = {
    1.43 +	1.e+36,
    1.44 +	1.e+36,
    1.45 +	1.e+36,
    1.46 +	1.e+36,
    1.47 +	1.e+36
    1.48 +};
    1.49 +static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars);
    1.50 +static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total ";
    1.51 +static char *time_name = "time", *lon_name = "lon", *lat_name = "lat";
    1.52 +static char *long_name_name = "long_name", *units_name = "units",
    1.53 +	*cell_method_name = "cell_method", *FillValue_name = "_FillValue",
    1.54 +	*missing_value_name = "missing_value";
    1.55 +static char *cell_method = "time: mean";
    1.56 +
    1.57 +void wrap_nc(int status)
    1.58 +{
    1.59 +	if (status != NC_NOERR) {
    1.60 +		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
    1.61 +		exit(1);
    1.62 +	}
    1.63 +
    1.64 +	return;
    1.65 +}
    1.66 +
    1.67 +void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1,
    1.68 +	size_t d2, char FillFlag, float FillValue, size_t tslice, float *var,
    1.69 +	double *tvar)
    1.70 +{
    1.71 +	int i, j;
    1.72 +	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
    1.73 +
    1.74 +	/* Initialize tvar to FillValue */
    1.75 +	for (i = 0; i < (d1 * d2); i++)
    1.76 +		tvar[i] = (double)FillValue;
    1.77 +
    1.78 +
    1.79 +	start[0] = tslice;
    1.80 +	count[0] = 1;
    1.81 +	start[2] = 0;
    1.82 +	count[2] = d1;
    1.83 +	start[3] = 0;
    1.84 +	count[3] = d2;
    1.85 +
    1.86 +	/* March through the levels, totalling as we go */
    1.87 +	for (j = 0; j < nlevels; j++) {
    1.88 +		start[1] = j;
    1.89 +		count[1] = 1;
    1.90 +		wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
    1.91 +		for (i = 0; i < (d1 * d2); i++) {
    1.92 +			if (tvar[i] == (double)FillValue)
    1.93 +				tvar[i] = (double)var[i];
    1.94 +			else
    1.95 +				tvar[i] += (double)var[i];
    1.96 +		}
    1.97 +	}
    1.98 +
    1.99 +	return;
   1.100 +}
   1.101 +
   1.102 +void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)
   1.103 +{
   1.104 +	int i;
   1.105 +	float **svars;
   1.106 +
   1.107 +	if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {
   1.108 +		perror("alloc_svars:svarsp");
   1.109 +		exit(2);
   1.110 +	}
   1.111 +	svars = *svarsp;
   1.112 +	for (i = 0; i < nsvars; i++) {
   1.113 +		if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) {
   1.114 +			perror("alloc_svars:svars[i]");
   1.115 +			exit(2);
   1.116 +		}
   1.117 +	}
   1.118 +	return;
   1.119 +}
   1.120 +
   1.121 +void free_svars(size_t nsvars, float **svars)
   1.122 +{
   1.123 +	int i;
   1.124 +
   1.125 +	for (i = 0; i < nsvars; i++)
   1.126 +		free(svars[i]);
   1.127 +	free(svars);
   1.128 +	return;
   1.129 +}
   1.130 +
   1.131 +void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,
   1.132 +	float *var)
   1.133 +{
   1.134 +	int varid;
   1.135 +	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
   1.136 +
   1.137 +	start[0] = tslice;
   1.138 +	count[0] = 1;
   1.139 +	start[1] = 0;
   1.140 +	count[1] = d1;
   1.141 +	start[2] = 0;
   1.142 +	count[2] = d2;
   1.143 +
   1.144 +	/* Get variable id */
   1.145 +	wrap_nc(nc_inq_varid(ncid, name, &varid));
   1.146 +	/* Assume the correct size and read the variable */
   1.147 +	wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
   1.148 +
   1.149 +	return;
   1.150 +}
   1.151 +
   1.152 +void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,
   1.153 +	int num, double *tvar)
   1.154 +{
   1.155 +	int i;
   1.156 +	size_t nsvars;
   1.157 +	float **svars; /* source variables */
   1.158 +	double denom;
   1.159 +
   1.160 +	/* Initialize tvar to FillValue */
   1.161 +	for (i = 0; i < (d1 * d2); i++)
   1.162 +		tvar[i] = (double)combo_FillValues[num];
   1.163 +
   1.164 +	switch (num) {
   1.165 +		case 0: /* LATENT */
   1.166 +			nsvars = 3;
   1.167 +			alloc_svars(nsvars, d1, d2, &svars);
   1.168 +			/* read timeslices */
   1.169 +			read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]);
   1.170 +			read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]);
   1.171 +			read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]);
   1.172 +			/* compute new variable */
   1.173 +			for (i = 0; i < (d1 * d2); i++) {
   1.174 +				if (svars[0][i] != combo_FillValues[num] &&
   1.175 +					svars[1][i] != combo_FillValues[num] &&
   1.176 +					svars[2][i] != combo_FillValues[num]) {
   1.177 +					tvar[i] = (double)svars[0][i]
   1.178 +						+ (double)svars[1][i]
   1.179 +						+ (double)svars[2][i];
   1.180 +				}
   1.181 +			}
   1.182 +			free_svars(nsvars, svars);
   1.183 +			break;
   1.184 +		case 1: /* ET */
   1.185 +			nsvars = 3;
   1.186 +			alloc_svars(nsvars, d1, d2, &svars);
   1.187 +			/* read timeslices */
   1.188 +			read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]);
   1.189 +			read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]);
   1.190 +			read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]);
   1.191 +			/* compute new variable */
   1.192 +			for (i = 0; i < (d1 * d2); i++) {
   1.193 +				if (svars[0][i] != combo_FillValues[num] &&
   1.194 +					svars[1][i] != combo_FillValues[num] &&
   1.195 +					svars[2][i] != combo_FillValues[num]) {
   1.196 +					tvar[i] = (double)svars[0][i]
   1.197 +						+ (double)svars[1][i]
   1.198 +						+ (double)svars[2][i];
   1.199 +				}
   1.200 +			}
   1.201 +			free_svars(nsvars, svars);
   1.202 +			break;
   1.203 +		case 2: /* ALL_SKY_ALBEDO */
   1.204 +			nsvars = 2;
   1.205 +			alloc_svars(nsvars, d1, d2, &svars);
   1.206 +			/* read timeslices */
   1.207 +			read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]);
   1.208 +			read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]);
   1.209 +			/* compute new variable */
   1.210 +			for (i = 0; i < (d1 * d2); i++) {
   1.211 +				if (svars[0][i] != combo_FillValues[num] &&
   1.212 +					svars[1][i] != combo_FillValues[num] &&
   1.213 +					(svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) {
   1.214 +					tvar[i] = (double)svars[0][i]
   1.215 +						/ (double)svars[1][i];
   1.216 +				}
   1.217 +			}
   1.218 +			free_svars(nsvars, svars);
   1.219 +			break;
   1.220 +		case 3: /* BLACK_SKY_ALBEDO */
   1.221 +			nsvars = 4;
   1.222 +			alloc_svars(nsvars, d1, d2, &svars);
   1.223 +			/* read timeslices */
   1.224 +			read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]);
   1.225 +			read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]);
   1.226 +			read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]);
   1.227 +			read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]);
   1.228 +			/* compute new variable */
   1.229 +			for (i = 0; i < (d1 * d2); i++) {
   1.230 +				if (svars[0][i] != combo_FillValues[num] &&
   1.231 +					svars[1][i] != combo_FillValues[num] &&
   1.232 +					svars[2][i] != combo_FillValues[num] &&
   1.233 +					svars[3][i] != combo_FillValues[num]) {
   1.234 +					denom = (double)svars[2][i] + (double)svars[3][i];
   1.235 +					if (denom > 1.e-35 || denom < -1.e-35) {
   1.236 +						tvar[i] = ((double)svars[0][i]
   1.237 +							+ (double)svars[1][i]) /
   1.238 +							((double)svars[2][i]
   1.239 +							+ (double)svars[3][i]);
   1.240 +					}
   1.241 +				}
   1.242 +			}
   1.243 +			free_svars(nsvars, svars);
   1.244 +			break;
   1.245 +		case 4: /* NETRAD */
   1.246 +			nsvars = 2;
   1.247 +			alloc_svars(nsvars, d1, d2, &svars);
   1.248 +			/* read timeslices */
   1.249 +			read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]);
   1.250 +			read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]);
   1.251 +			/* compute new variable */
   1.252 +			for (i = 0; i < (d1 * d2); i++) {
   1.253 +				if (svars[0][i] != combo_FillValues[num] &&
   1.254 +					svars[1][i] != combo_FillValues[num]) {
   1.255 +					tvar[i] = (double)svars[0][i]
   1.256 +						- (double)svars[1][i];
   1.257 +				}
   1.258 +			}
   1.259 +			free_svars(nsvars, svars);
   1.260 +			break;
   1.261 +		default:
   1.262 +			fprintf(stderr, "Unknown combination variable %d!\n", num);
   1.263 +			exit(3);
   1.264 +	}
   1.265 +	
   1.266 +	return;
   1.267 +}
   1.268 +
   1.269 +void write_timeslice(int ncid, int varid, size_t d1, size_t d2,
   1.270 +	size_t tslice, float *var)
   1.271 +{
   1.272 +	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
   1.273 +
   1.274 +	start[0] = tslice;
   1.275 +	count[0] = 1;
   1.276 +	start[1] = 0;
   1.277 +	count[1] = d1;
   1.278 +	start[2] = 0;
   1.279 +	count[2] = d2;
   1.280 +
   1.281 +	wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));
   1.282 +
   1.283 +	return;
   1.284 +}
   1.285 +
   1.286 +void usage(char *program)
   1.287 +{
   1.288 +	fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",
   1.289 +		program);
   1.290 +	fprintf(stderr, "Requires at least one history output file name\n");
   1.291 +	exit(5);
   1.292 +}
   1.293 +
   1.294 +int main(int argc, char **argv)
   1.295 +{
   1.296 +	char **ifnames;
   1.297 +	int i, j, k, nifnames;
   1.298 +	int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts,
   1.299 +		dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS],
   1.300 +		new_varid, lat_dimid, lon_dimid;
   1.301 +	nc_type xtype, atype;
   1.302 +	char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1],
   1.303 +		att_name[NC_MAX_NAME+1], *longname, *new_longname;
   1.304 +	size_t tlen, ts, alen, nlevels, d1, d2;
   1.305 +	float FillValue;
   1.306 +	char FillFlag;
   1.307 +	float *var;
   1.308 +	double *tvar;
   1.309 +
   1.310 +#ifdef DEBUG
   1.311 +	printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);
   1.312 +	printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);
   1.313 +#endif /* DEBUG */
   1.314 +
   1.315 +	/* Require at least one parameter (the file on which to work) */
   1.316 +	if (argc < 2)
   1.317 +		usage(argv[0]);
   1.318 +
   1.319 +	/* Put the list of filenames in an array of strings */
   1.320 +	if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {
   1.321 +		perror("ifnames");
   1.322 +		exit(2);
   1.323 +	}
   1.324 +	nifnames = 0;
   1.325 +	for (i = 1; i < argc; i++) {
   1.326 +		if (!(ifnames[nifnames++] = strdup(argv[i]))) {
   1.327 +			perror("ifnames[nifnames++]");
   1.328 +			exit(2);
   1.329 +		}
   1.330 +	}
   1.331 +
   1.332 +	/*
   1.333 +	 * March through all input files adding new variables for totals
   1.334 +	 * across levels.
   1.335 +	 */
   1.336 +	for (i = 0; i < nifnames; i++) {
   1.337 +		printf("Processing %s\n", ifnames[i]);
   1.338 +		/* Open file */
   1.339 +		wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));
   1.340 +		/*
   1.341 +		 * Inquire about number of dimensions, variables, global
   1.342 +		 * attributes and the ID of the unlimited dimension
   1.343 +		 */
   1.344 +		wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid));
   1.345 +		wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen));
   1.346 +		if (strcmp(name, time_name)) {
   1.347 +			fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]);
   1.348 +			exit(4);
   1.349 +		}
   1.350 +		/* Retrieve lat and lon dimension IDs */
   1.351 +		wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid));
   1.352 +		wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid));
   1.353 +		/* Put file into define mode */
   1.354 +		wrap_nc(nc_redef(ncid));
   1.355 +		/*
   1.356 +		 * Define all of the new variables first to improve
   1.357 +		 * performance.  This means some calls are made twice.
   1.358 +		 */
   1.359 +		/* First, add the new total variables */
   1.360 +		for (j = 0; j < nmtotal_vars; j++) {
   1.361 +			/* Figure out source variable information */
   1.362 +			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
   1.363 +			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
   1.364 +			/* Make sure it has the correct number of dimensions */
   1.365 +			if (ndims != 4) {
   1.366 +				fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);
   1.367 +				exit(4);
   1.368 +			}
   1.369 +			/* Make sure time is the first dimension */
   1.370 +			if (dimids[0] != unlimdimid) {
   1.371 +				fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name);
   1.372 +				exit(4);
   1.373 +			}
   1.374 +			/* Make sure it is a float type */
   1.375 +			if (xtype != NC_FLOAT) {
   1.376 +				fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");
   1.377 +				exit(4);
   1.378 +			}
   1.379 +			/* Set dimensions for new variable */
   1.380 +			new_ndims = 3;
   1.381 +			new_dimids[0] = unlimdimid;
   1.382 +			new_dimids[1] = dimids[2];
   1.383 +			new_dimids[2] = dimids[3];
   1.384 +			/* Define new variable */
   1.385 +			sprintf(new_name, "%s%s", varname_total_prefix, name);
   1.386 +			printf("\tAdding variable %s\n", new_name);
   1.387 +			wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims,
   1.388 +				new_dimids, &new_varid));
   1.389 +			/* Copy attributes from source variable to new one */
   1.390 +			for (k = 0; k < natts; k++) {
   1.391 +				wrap_nc(nc_inq_attname(ncid, varid, k, att_name));
   1.392 +				if (!strcmp(att_name, "long_name")) {
   1.393 +					wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen));
   1.394 +					if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) {
   1.395 +						perror("longname");
   1.396 +						exit(2);
   1.397 +					}
   1.398 +					wrap_nc(nc_get_att_text(ncid, varid, att_name, longname));
   1.399 +					longname[alen] = '\0';
   1.400 +					if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) {
   1.401 +						perror("new_longname");
   1.402 +						exit(2);
   1.403 +					}
   1.404 +					sprintf(new_longname, "%s%s", longname_total_prefix, longname);
   1.405 +					wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname));
   1.406 +					free(new_longname);
   1.407 +					free(longname);
   1.408 +				}
   1.409 +				else
   1.410 +					wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));
   1.411 +			}
   1.412 +		}
   1.413 +		/* Second, add the new total variables */
   1.414 +		for (j = 0; j < nmcombo_vars; j++) {
   1.415 +			/* Set dimensions for new variable */
   1.416 +			new_ndims = 3;
   1.417 +			new_dimids[0] = unlimdimid;
   1.418 +			new_dimids[1] = lat_dimid;
   1.419 +			new_dimids[2] = lon_dimid;
   1.420 +			/* Define new variable */
   1.421 +			printf("\tAdding variable %s\n", combo_vars[j]);
   1.422 +			wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT,
   1.423 +				new_ndims, new_dimids, &new_varid));
   1.424 +			/* Add attributes */
   1.425 +			wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name,
   1.426 +				strlen(combo_long_names[j]),
   1.427 +				combo_long_names[j]));
   1.428 +			wrap_nc(nc_put_att_text(ncid, new_varid, units_name,
   1.429 +				strlen(combo_units[j]), combo_units[j]));
   1.430 +			wrap_nc(nc_put_att_text(ncid, new_varid,
   1.431 +				cell_method_name, strlen(cell_method),
   1.432 +				cell_method));
   1.433 +			wrap_nc(nc_put_att_float(ncid, new_varid,
   1.434 +				FillValue_name, NC_FLOAT, 1,
   1.435 +				&combo_FillValues[j]));
   1.436 +			wrap_nc(nc_put_att_float(ncid, new_varid,
   1.437 +				missing_value_name, NC_FLOAT, 1,
   1.438 +				&combo_FillValues[j]));
   1.439 +		}
   1.440 +		/* Leave define mode */
   1.441 +		wrap_nc(nc_enddef(ncid));
   1.442 +		/* Now actually compute and write out the new total variables */
   1.443 +		for (j = 0; j < nmtotal_vars; j++) {
   1.444 +			/* Figure out source variable information */
   1.445 +			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
   1.446 +			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
   1.447 +			/* Check for _FillValue */
   1.448 +			FillFlag = 0;
   1.449 +			FillValue = 0.;
   1.450 +			if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) {
   1.451 +				if (atype == NC_FLOAT && alen) {
   1.452 +					wrap_nc(nc_get_att_float(ncid, varid,
   1.453 +						FillValue_name, &FillValue));
   1.454 +					FillFlag = 1;
   1.455 +				}
   1.456 +			}
   1.457 +			/* Set dimensions for new variable */
   1.458 +			new_ndims = 3;
   1.459 +			new_dimids[0] = unlimdimid;
   1.460 +			new_dimids[1] = dimids[2];
   1.461 +			new_dimids[2] = dimids[3];
   1.462 +
   1.463 +			sprintf(new_name, "%s%s", varname_total_prefix, name);
   1.464 +			printf("\tComputing and writing %s\n", new_name);
   1.465 +			/* Retrieve the new varid again */
   1.466 +			wrap_nc(nc_inq_varid(ncid, new_name, &new_varid));
   1.467 +			/*
   1.468 +			 * Figure out dimensions and total size
   1.469 +			 * for a single time slice
   1.470 +			 */
   1.471 +			wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels));
   1.472 +			wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1));
   1.473 +			wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2));
   1.474 +
   1.475 +			/*
   1.476 +			 * Allocate space for reading in time slice and for
   1.477 +			 * the sum
   1.478 +			 */
   1.479 +			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
   1.480 +				perror("var");
   1.481 +				exit(2);
   1.482 +			}
   1.483 +			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
   1.484 +				perror("tvar");
   1.485 +				exit(2);
   1.486 +			}
   1.487 +
   1.488 +			/* Read timeslices, total, and write out new timeslices */
   1.489 +			for (ts = 0; ts < tlen; ts++) {
   1.490 +				read_total_timeslice(ncid, varid,
   1.491 +					nlevels, d1, d2, FillFlag, FillValue,
   1.492 +					ts, var, tvar);
   1.493 +				for (k = 0; k < (d1*d2); k++)
   1.494 +					var[k] = (float)tvar[k];
   1.495 +				write_timeslice(ncid, new_varid, d1, d2,
   1.496 +					ts, var);
   1.497 +			}
   1.498 +
   1.499 +			free(var);
   1.500 +			free(tvar);
   1.501 +		}
   1.502 +		/* Now actually compute and write out the new combo variables */
   1.503 +		for (j = 0; j < nmcombo_vars; j++) {
   1.504 +			/* Set dimensions for new variable */
   1.505 +			new_ndims = 3;
   1.506 +			new_dimids[0] = unlimdimid;
   1.507 +			new_dimids[1] = lat_dimid;
   1.508 +			new_dimids[2] = lon_dimid;
   1.509 +
   1.510 +			printf("\tComputing and writing %s\n", combo_vars[j]);
   1.511 +			/* Retrieve the new varid again */
   1.512 +			wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid));
   1.513 +			/*
   1.514 +			 * Figure out dimensions and total size
   1.515 +			 * for a single time slice
   1.516 +			 */
   1.517 +			wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));
   1.518 +			wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));
   1.519 +
   1.520 +			/*
   1.521 +			 * Allocate space for reading in time slice and for
   1.522 +			 * the sum
   1.523 +			 */
   1.524 +			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
   1.525 +				perror("var");
   1.526 +				exit(2);
   1.527 +			}
   1.528 +			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
   1.529 +				perror("tvar");
   1.530 +				exit(2);
   1.531 +			}
   1.532 +
   1.533 +			/* Read timeslices, compute new variable, and write */
   1.534 +			for (ts = 0; ts < tlen; ts++) {
   1.535 +				read_compute_timeslice(ncid, d1, d2, ts,
   1.536 +					j, tvar);
   1.537 +				for (k = 0; k < (d1*d2); k++)
   1.538 +					var[k] = (float)tvar[k];
   1.539 +				write_timeslice(ncid, new_varid, d1, d2,
   1.540 +					ts, var);
   1.541 +			}
   1.542 +
   1.543 +			free(var);
   1.544 +			free(tvar);
   1.545 +		}
   1.546 +		/* Close file */
   1.547 +		wrap_nc(nc_close(ncid));
   1.548 +	}
   1.549 +
   1.550 +	return 0;
   1.551 +}