forrest@0: #include <stdio.h>
forrest@0: #include <stdlib.h>
forrest@0: #include <string.h>
forrest@0: #include <math.h>
forrest@0: #include "netcdf.h"
forrest@0: 
forrest@0: /*
forrest@0:  * Loop through all history tapes from the Community Land Model adding new
forrest@0:  * fields that are totals of the multi-level fields included in total_vars[]
forrest@0:  * below.
forrest@0:  */
forrest@0: 
forrest@0: static char *total_vars[] = {
forrest@0: 	"SOILICE",
forrest@0: 	"SOILLIQ"
forrest@0: };
forrest@0: static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars);
forrest@0: static char *combo_vars[] = {
forrest@0: 	"LATENT",		/* FCTR + FCEV + FGEV */
forrest@0: 	"ET",			/* QVEGE + QVEGT + QSOIL */
forrest@0: 	"ALL_SKY_ALBEDO",	/* FSR / FSDS */
forrest@0: 	"BLACK_SKY_ALBEDO",	/* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */
forrest@0: 	"NETRAD",		/* FSA - FIRA */
forrest@0: };
forrest@0: static char *combo_long_names[] = {
forrest@0: 	"latent heat flux",
forrest@0: 	"evapotranspiration",
forrest@0: 	"surface all-sky albedo",
forrest@0: 	"surface black-sky albedo",
forrest@0: 	"net radiation"
forrest@0: };
forrest@0: static char *combo_units[] = {
forrest@0: 	"watt/m^2",
forrest@0: 	"m/s",
forrest@0: 	"unitless",
forrest@0: 	"unitless",
forrest@0: 	"watt/m^2"
forrest@0: };
forrest@0: static float combo_FillValues[] = {
forrest@0: 	1.e+36,
forrest@0: 	1.e+36,
forrest@0: 	1.e+36,
forrest@0: 	1.e+36,
forrest@0: 	1.e+36
forrest@0: };
forrest@0: static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars);
forrest@0: static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total ";
forrest@0: static char *time_name = "time", *lon_name = "lon", *lat_name = "lat";
forrest@0: static char *long_name_name = "long_name", *units_name = "units",
forrest@0: 	*cell_method_name = "cell_method", *FillValue_name = "_FillValue",
forrest@0: 	*missing_value_name = "missing_value";
forrest@0: static char *cell_method = "time: mean";
forrest@0: 
forrest@0: void wrap_nc(int status)
forrest@0: {
forrest@0: 	if (status != NC_NOERR) {
forrest@0: 		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
forrest@0: 		exit(1);
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1,
forrest@0: 	size_t d2, char FillFlag, float FillValue, size_t tslice, float *var,
forrest@0: 	double *tvar)
forrest@0: {
forrest@0: 	int i, j;
forrest@0: 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 
forrest@0: 	/* Initialize tvar to FillValue */
forrest@0: 	for (i = 0; i < (d1 * d2); i++)
forrest@0: 		tvar[i] = (double)FillValue;
forrest@0: 
forrest@0: 
forrest@0: 	start[0] = tslice;
forrest@0: 	count[0] = 1;
forrest@0: 	start[2] = 0;
forrest@0: 	count[2] = d1;
forrest@0: 	start[3] = 0;
forrest@0: 	count[3] = d2;
forrest@0: 
forrest@0: 	/* March through the levels, totalling as we go */
forrest@0: 	for (j = 0; j < nlevels; j++) {
forrest@0: 		start[1] = j;
forrest@0: 		count[1] = 1;
forrest@0: 		wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
forrest@0: 		for (i = 0; i < (d1 * d2); i++) {
forrest@0: 			if (tvar[i] == (double)FillValue)
forrest@0: 				tvar[i] = (double)var[i];
forrest@0: 			else
forrest@0: 				tvar[i] += (double)var[i];
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	float **svars;
forrest@0: 
forrest@0: 	if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {
forrest@0: 		perror("alloc_svars:svarsp");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	svars = *svarsp;
forrest@0: 	for (i = 0; i < nsvars; i++) {
forrest@0: 		if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) {
forrest@0: 			perror("alloc_svars:svars[i]");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 	}
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void free_svars(size_t nsvars, float **svars)
forrest@0: {
forrest@0: 	int i;
forrest@0: 
forrest@0: 	for (i = 0; i < nsvars; i++)
forrest@0: 		free(svars[i]);
forrest@0: 	free(svars);
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,
forrest@0: 	float *var)
forrest@0: {
forrest@0: 	int varid;
forrest@0: 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 
forrest@0: 	start[0] = tslice;
forrest@0: 	count[0] = 1;
forrest@0: 	start[1] = 0;
forrest@0: 	count[1] = d1;
forrest@0: 	start[2] = 0;
forrest@0: 	count[2] = d2;
forrest@0: 
forrest@0: 	/* Get variable id */
forrest@0: 	wrap_nc(nc_inq_varid(ncid, name, &varid));
forrest@0: 	/* Assume the correct size and read the variable */
forrest@0: 	wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,
forrest@0: 	int num, double *tvar)
forrest@0: {
forrest@0: 	int i;
forrest@0: 	size_t nsvars;
forrest@0: 	float **svars; /* source variables */
forrest@0: 	double denom;
forrest@0: 
forrest@0: 	/* Initialize tvar to FillValue */
forrest@0: 	for (i = 0; i < (d1 * d2); i++)
forrest@0: 		tvar[i] = (double)combo_FillValues[num];
forrest@0: 
forrest@0: 	switch (num) {
forrest@0: 		case 0: /* LATENT */
forrest@0: 			nsvars = 3;
forrest@0: 			alloc_svars(nsvars, d1, d2, &svars);
forrest@0: 			/* read timeslices */
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]);
forrest@0: 			/* compute new variable */
forrest@0: 			for (i = 0; i < (d1 * d2); i++) {
forrest@0: 				if (svars[0][i] != combo_FillValues[num] &&
forrest@0: 					svars[1][i] != combo_FillValues[num] &&
forrest@0: 					svars[2][i] != combo_FillValues[num]) {
forrest@0: 					tvar[i] = (double)svars[0][i]
forrest@0: 						+ (double)svars[1][i]
forrest@0: 						+ (double)svars[2][i];
forrest@0: 				}
forrest@0: 			}
forrest@0: 			free_svars(nsvars, svars);
forrest@0: 			break;
forrest@0: 		case 1: /* ET */
forrest@0: 			nsvars = 3;
forrest@0: 			alloc_svars(nsvars, d1, d2, &svars);
forrest@0: 			/* read timeslices */
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]);
forrest@0: 			/* compute new variable */
forrest@0: 			for (i = 0; i < (d1 * d2); i++) {
forrest@0: 				if (svars[0][i] != combo_FillValues[num] &&
forrest@0: 					svars[1][i] != combo_FillValues[num] &&
forrest@0: 					svars[2][i] != combo_FillValues[num]) {
forrest@0: 					tvar[i] = (double)svars[0][i]
forrest@0: 						+ (double)svars[1][i]
forrest@0: 						+ (double)svars[2][i];
forrest@0: 				}
forrest@0: 			}
forrest@0: 			free_svars(nsvars, svars);
forrest@0: 			break;
forrest@0: 		case 2: /* ALL_SKY_ALBEDO */
forrest@0: 			nsvars = 2;
forrest@0: 			alloc_svars(nsvars, d1, d2, &svars);
forrest@0: 			/* read timeslices */
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]);
forrest@0: 			/* compute new variable */
forrest@0: 			for (i = 0; i < (d1 * d2); i++) {
forrest@0: 				if (svars[0][i] != combo_FillValues[num] &&
forrest@0: 					svars[1][i] != combo_FillValues[num] &&
forrest@0: 					(svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) {
forrest@0: 					tvar[i] = (double)svars[0][i]
forrest@0: 						/ (double)svars[1][i];
forrest@0: 				}
forrest@0: 			}
forrest@0: 			free_svars(nsvars, svars);
forrest@0: 			break;
forrest@0: 		case 3: /* BLACK_SKY_ALBEDO */
forrest@0: 			nsvars = 4;
forrest@0: 			alloc_svars(nsvars, d1, d2, &svars);
forrest@0: 			/* read timeslices */
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]);
forrest@0: 			/* compute new variable */
forrest@0: 			for (i = 0; i < (d1 * d2); i++) {
forrest@0: 				if (svars[0][i] != combo_FillValues[num] &&
forrest@0: 					svars[1][i] != combo_FillValues[num] &&
forrest@0: 					svars[2][i] != combo_FillValues[num] &&
forrest@0: 					svars[3][i] != combo_FillValues[num]) {
forrest@0: 					denom = (double)svars[2][i] + (double)svars[3][i];
forrest@0: 					if (denom > 1.e-35 || denom < -1.e-35) {
forrest@0: 						tvar[i] = ((double)svars[0][i]
forrest@0: 							+ (double)svars[1][i]) /
forrest@0: 							((double)svars[2][i]
forrest@0: 							+ (double)svars[3][i]);
forrest@0: 					}
forrest@0: 				}
forrest@0: 			}
forrest@0: 			free_svars(nsvars, svars);
forrest@0: 			break;
forrest@0: 		case 4: /* NETRAD */
forrest@0: 			nsvars = 2;
forrest@0: 			alloc_svars(nsvars, d1, d2, &svars);
forrest@0: 			/* read timeslices */
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]);
forrest@0: 			read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]);
forrest@0: 			/* compute new variable */
forrest@0: 			for (i = 0; i < (d1 * d2); i++) {
forrest@0: 				if (svars[0][i] != combo_FillValues[num] &&
forrest@0: 					svars[1][i] != combo_FillValues[num]) {
forrest@0: 					tvar[i] = (double)svars[0][i]
forrest@0: 						- (double)svars[1][i];
forrest@0: 				}
forrest@0: 			}
forrest@0: 			free_svars(nsvars, svars);
forrest@0: 			break;
forrest@0: 		default:
forrest@0: 			fprintf(stderr, "Unknown combination variable %d!\n", num);
forrest@0: 			exit(3);
forrest@0: 	}
forrest@0: 	
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void write_timeslice(int ncid, int varid, size_t d1, size_t d2,
forrest@0: 	size_t tslice, float *var)
forrest@0: {
forrest@0: 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
forrest@0: 
forrest@0: 	start[0] = tslice;
forrest@0: 	count[0] = 1;
forrest@0: 	start[1] = 0;
forrest@0: 	count[1] = d1;
forrest@0: 	start[2] = 0;
forrest@0: 	count[2] = d2;
forrest@0: 
forrest@0: 	wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));
forrest@0: 
forrest@0: 	return;
forrest@0: }
forrest@0: 
forrest@0: void usage(char *program)
forrest@0: {
forrest@0: 	fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",
forrest@0: 		program);
forrest@0: 	fprintf(stderr, "Requires at least one history output file name\n");
forrest@0: 	exit(5);
forrest@0: }
forrest@0: 
forrest@0: int main(int argc, char **argv)
forrest@0: {
forrest@0: 	char **ifnames;
forrest@0: 	int i, j, k, nifnames;
forrest@0: 	int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts,
forrest@0: 		dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS],
forrest@0: 		new_varid, lat_dimid, lon_dimid;
forrest@0: 	nc_type xtype, atype;
forrest@0: 	char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1],
forrest@0: 		att_name[NC_MAX_NAME+1], *longname, *new_longname;
forrest@0: 	size_t tlen, ts, alen, nlevels, d1, d2;
forrest@0: 	float FillValue;
forrest@0: 	char FillFlag;
forrest@0: 	float *var;
forrest@0: 	double *tvar;
forrest@0: 
forrest@0: #ifdef DEBUG
forrest@0: 	printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);
forrest@0: 	printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);
forrest@0: #endif /* DEBUG */
forrest@0: 
forrest@0: 	/* Require at least one parameter (the file on which to work) */
forrest@0: 	if (argc < 2)
forrest@0: 		usage(argv[0]);
forrest@0: 
forrest@0: 	/* Put the list of filenames in an array of strings */
forrest@0: 	if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {
forrest@0: 		perror("ifnames");
forrest@0: 		exit(2);
forrest@0: 	}
forrest@0: 	nifnames = 0;
forrest@0: 	for (i = 1; i < argc; i++) {
forrest@0: 		if (!(ifnames[nifnames++] = strdup(argv[i]))) {
forrest@0: 			perror("ifnames[nifnames++]");
forrest@0: 			exit(2);
forrest@0: 		}
forrest@0: 	}
forrest@0: 
forrest@0: 	/*
forrest@0: 	 * March through all input files adding new variables for totals
forrest@0: 	 * across levels.
forrest@0: 	 */
forrest@0: 	for (i = 0; i < nifnames; i++) {
forrest@0: 		printf("Processing %s\n", ifnames[i]);
forrest@0: 		/* Open file */
forrest@0: 		wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));
forrest@0: 		/*
forrest@0: 		 * Inquire about number of dimensions, variables, global
forrest@0: 		 * attributes and the ID of the unlimited dimension
forrest@0: 		 */
forrest@0: 		wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid));
forrest@0: 		wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen));
forrest@0: 		if (strcmp(name, time_name)) {
forrest@0: 			fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]);
forrest@0: 			exit(4);
forrest@0: 		}
forrest@0: 		/* Retrieve lat and lon dimension IDs */
forrest@0: 		wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid));
forrest@0: 		wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid));
forrest@0: 		/* Put file into define mode */
forrest@0: 		wrap_nc(nc_redef(ncid));
forrest@0: 		/*
forrest@0: 		 * Define all of the new variables first to improve
forrest@0: 		 * performance.  This means some calls are made twice.
forrest@0: 		 */
forrest@0: 		/* First, add the new total variables */
forrest@0: 		for (j = 0; j < nmtotal_vars; j++) {
forrest@0: 			/* Figure out source variable information */
forrest@0: 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
forrest@0: 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
forrest@0: 			/* Make sure it has the correct number of dimensions */
forrest@0: 			if (ndims != 4) {
forrest@0: 				fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);
forrest@0: 				exit(4);
forrest@0: 			}
forrest@0: 			/* Make sure time is the first dimension */
forrest@0: 			if (dimids[0] != unlimdimid) {
forrest@0: 				fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name);
forrest@0: 				exit(4);
forrest@0: 			}
forrest@0: 			/* Make sure it is a float type */
forrest@0: 			if (xtype != NC_FLOAT) {
forrest@0: 				fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");
forrest@0: 				exit(4);
forrest@0: 			}
forrest@0: 			/* Set dimensions for new variable */
forrest@0: 			new_ndims = 3;
forrest@0: 			new_dimids[0] = unlimdimid;
forrest@0: 			new_dimids[1] = dimids[2];
forrest@0: 			new_dimids[2] = dimids[3];
forrest@0: 			/* Define new variable */
forrest@0: 			sprintf(new_name, "%s%s", varname_total_prefix, name);
forrest@0: 			printf("\tAdding variable %s\n", new_name);
forrest@0: 			wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims,
forrest@0: 				new_dimids, &new_varid));
forrest@0: 			/* Copy attributes from source variable to new one */
forrest@0: 			for (k = 0; k < natts; k++) {
forrest@0: 				wrap_nc(nc_inq_attname(ncid, varid, k, att_name));
forrest@0: 				if (!strcmp(att_name, "long_name")) {
forrest@0: 					wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen));
forrest@0: 					if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) {
forrest@0: 						perror("longname");
forrest@0: 						exit(2);
forrest@0: 					}
forrest@0: 					wrap_nc(nc_get_att_text(ncid, varid, att_name, longname));
forrest@0: 					longname[alen] = '\0';
forrest@0: 					if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) {
forrest@0: 						perror("new_longname");
forrest@0: 						exit(2);
forrest@0: 					}
forrest@0: 					sprintf(new_longname, "%s%s", longname_total_prefix, longname);
forrest@0: 					wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname));
forrest@0: 					free(new_longname);
forrest@0: 					free(longname);
forrest@0: 				}
forrest@0: 				else
forrest@0: 					wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));
forrest@0: 			}
forrest@0: 		}
forrest@0: 		/* Second, add the new total variables */
forrest@0: 		for (j = 0; j < nmcombo_vars; j++) {
forrest@0: 			/* Set dimensions for new variable */
forrest@0: 			new_ndims = 3;
forrest@0: 			new_dimids[0] = unlimdimid;
forrest@0: 			new_dimids[1] = lat_dimid;
forrest@0: 			new_dimids[2] = lon_dimid;
forrest@0: 			/* Define new variable */
forrest@0: 			printf("\tAdding variable %s\n", combo_vars[j]);
forrest@0: 			wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT,
forrest@0: 				new_ndims, new_dimids, &new_varid));
forrest@0: 			/* Add attributes */
forrest@0: 			wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name,
forrest@0: 				strlen(combo_long_names[j]),
forrest@0: 				combo_long_names[j]));
forrest@0: 			wrap_nc(nc_put_att_text(ncid, new_varid, units_name,
forrest@0: 				strlen(combo_units[j]), combo_units[j]));
forrest@0: 			wrap_nc(nc_put_att_text(ncid, new_varid,
forrest@0: 				cell_method_name, strlen(cell_method),
forrest@0: 				cell_method));
forrest@0: 			wrap_nc(nc_put_att_float(ncid, new_varid,
forrest@0: 				FillValue_name, NC_FLOAT, 1,
forrest@0: 				&combo_FillValues[j]));
forrest@0: 			wrap_nc(nc_put_att_float(ncid, new_varid,
forrest@0: 				missing_value_name, NC_FLOAT, 1,
forrest@0: 				&combo_FillValues[j]));
forrest@0: 		}
forrest@0: 		/* Leave define mode */
forrest@0: 		wrap_nc(nc_enddef(ncid));
forrest@0: 		/* Now actually compute and write out the new total variables */
forrest@0: 		for (j = 0; j < nmtotal_vars; j++) {
forrest@0: 			/* Figure out source variable information */
forrest@0: 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
forrest@0: 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
forrest@0: 			/* Check for _FillValue */
forrest@0: 			FillFlag = 0;
forrest@0: 			FillValue = 0.;
forrest@0: 			if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) {
forrest@0: 				if (atype == NC_FLOAT && alen) {
forrest@0: 					wrap_nc(nc_get_att_float(ncid, varid,
forrest@0: 						FillValue_name, &FillValue));
forrest@0: 					FillFlag = 1;
forrest@0: 				}
forrest@0: 			}
forrest@0: 			/* Set dimensions for new variable */
forrest@0: 			new_ndims = 3;
forrest@0: 			new_dimids[0] = unlimdimid;
forrest@0: 			new_dimids[1] = dimids[2];
forrest@0: 			new_dimids[2] = dimids[3];
forrest@0: 
forrest@0: 			sprintf(new_name, "%s%s", varname_total_prefix, name);
forrest@0: 			printf("\tComputing and writing %s\n", new_name);
forrest@0: 			/* Retrieve the new varid again */
forrest@0: 			wrap_nc(nc_inq_varid(ncid, new_name, &new_varid));
forrest@0: 			/*
forrest@0: 			 * Figure out dimensions and total size
forrest@0: 			 * for a single time slice
forrest@0: 			 */
forrest@0: 			wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels));
forrest@0: 			wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1));
forrest@0: 			wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2));
forrest@0: 
forrest@0: 			/*
forrest@0: 			 * Allocate space for reading in time slice and for
forrest@0: 			 * the sum
forrest@0: 			 */
forrest@0: 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
forrest@0: 				perror("var");
forrest@0: 				exit(2);
forrest@0: 			}
forrest@0: 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
forrest@0: 				perror("tvar");
forrest@0: 				exit(2);
forrest@0: 			}
forrest@0: 
forrest@0: 			/* Read timeslices, total, and write out new timeslices */
forrest@0: 			for (ts = 0; ts < tlen; ts++) {
forrest@0: 				read_total_timeslice(ncid, varid,
forrest@0: 					nlevels, d1, d2, FillFlag, FillValue,
forrest@0: 					ts, var, tvar);
forrest@0: 				for (k = 0; k < (d1*d2); k++)
forrest@0: 					var[k] = (float)tvar[k];
forrest@0: 				write_timeslice(ncid, new_varid, d1, d2,
forrest@0: 					ts, var);
forrest@0: 			}
forrest@0: 
forrest@0: 			free(var);
forrest@0: 			free(tvar);
forrest@0: 		}
forrest@0: 		/* Now actually compute and write out the new combo variables */
forrest@0: 		for (j = 0; j < nmcombo_vars; j++) {
forrest@0: 			/* Set dimensions for new variable */
forrest@0: 			new_ndims = 3;
forrest@0: 			new_dimids[0] = unlimdimid;
forrest@0: 			new_dimids[1] = lat_dimid;
forrest@0: 			new_dimids[2] = lon_dimid;
forrest@0: 
forrest@0: 			printf("\tComputing and writing %s\n", combo_vars[j]);
forrest@0: 			/* Retrieve the new varid again */
forrest@0: 			wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid));
forrest@0: 			/*
forrest@0: 			 * Figure out dimensions and total size
forrest@0: 			 * for a single time slice
forrest@0: 			 */
forrest@0: 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));
forrest@0: 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));
forrest@0: 
forrest@0: 			/*
forrest@0: 			 * Allocate space for reading in time slice and for
forrest@0: 			 * the sum
forrest@0: 			 */
forrest@0: 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
forrest@0: 				perror("var");
forrest@0: 				exit(2);
forrest@0: 			}
forrest@0: 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
forrest@0: 				perror("tvar");
forrest@0: 				exit(2);
forrest@0: 			}
forrest@0: 
forrest@0: 			/* Read timeslices, compute new variable, and write */
forrest@0: 			for (ts = 0; ts < tlen; ts++) {
forrest@0: 				read_compute_timeslice(ncid, d1, d2, ts,
forrest@0: 					j, tvar);
forrest@0: 				for (k = 0; k < (d1*d2); k++)
forrest@0: 					var[k] = (float)tvar[k];
forrest@0: 				write_timeslice(ncid, new_varid, d1, d2,
forrest@0: 					ts, var);
forrest@0: 			}
forrest@0: 
forrest@0: 			free(var);
forrest@0: 			free(tvar);
forrest@0: 		}
forrest@0: 		/* Close file */
forrest@0: 		wrap_nc(nc_close(ncid));
forrest@0: 	}
forrest@0: 
forrest@0: 	return 0;
forrest@0: }