add_total_fields.c
changeset 1 2ce4ee911439
equal deleted inserted replaced
-1:000000000000 0:c9c5bb64e022
       
     1 #include <stdio.h>
       
     2 #include <stdlib.h>
       
     3 #include <string.h>
       
     4 #include <math.h>
       
     5 #include "netcdf.h"
       
     6 
       
     7 /*
       
     8  * Loop through all history tapes from the Community Land Model adding new
       
     9  * fields that are totals of the multi-level fields included in total_vars[]
       
    10  * below.
       
    11  */
       
    12 
       
    13 static char *total_vars[] = {
       
    14 	"SOILICE",
       
    15 	"SOILLIQ"
       
    16 };
       
    17 static int nmtotal_vars = sizeof(total_vars)/sizeof(*total_vars);
       
    18 static char *combo_vars[] = {
       
    19 	"LATENT",		/* FCTR + FCEV + FGEV */
       
    20 	"ET",			/* QVEGE + QVEGT + QSOIL */
       
    21 	"ALL_SKY_ALBEDO",	/* FSR / FSDS */
       
    22 	"BLACK_SKY_ALBEDO",	/* (FSRNDLN + FSRVDLN)/(FSDSNDLN + FSDSVDLN) */
       
    23 	"NETRAD",		/* FSA - FIRA */
       
    24 };
       
    25 static char *combo_long_names[] = {
       
    26 	"latent heat flux",
       
    27 	"evapotranspiration",
       
    28 	"surface all-sky albedo",
       
    29 	"surface black-sky albedo",
       
    30 	"net radiation"
       
    31 };
       
    32 static char *combo_units[] = {
       
    33 	"watt/m^2",
       
    34 	"m/s",
       
    35 	"unitless",
       
    36 	"unitless",
       
    37 	"watt/m^2"
       
    38 };
       
    39 static float combo_FillValues[] = {
       
    40 	1.e+36,
       
    41 	1.e+36,
       
    42 	1.e+36,
       
    43 	1.e+36,
       
    44 	1.e+36
       
    45 };
       
    46 static int nmcombo_vars = sizeof(combo_vars)/sizeof(*combo_vars);
       
    47 static char *varname_total_prefix = "TOTAL_", *longname_total_prefix = "total ";
       
    48 static char *time_name = "time", *lon_name = "lon", *lat_name = "lat";
       
    49 static char *long_name_name = "long_name", *units_name = "units",
       
    50 	*cell_method_name = "cell_method", *FillValue_name = "_FillValue",
       
    51 	*missing_value_name = "missing_value";
       
    52 static char *cell_method = "time: mean";
       
    53 
       
    54 void wrap_nc(int status)
       
    55 {
       
    56 	if (status != NC_NOERR) {
       
    57 		fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
       
    58 		exit(1);
       
    59 	}
       
    60 
       
    61 	return;
       
    62 }
       
    63 
       
    64 void read_total_timeslice(int ncid, int varid, size_t nlevels, size_t d1,
       
    65 	size_t d2, char FillFlag, float FillValue, size_t tslice, float *var,
       
    66 	double *tvar)
       
    67 {
       
    68 	int i, j;
       
    69 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
       
    70 
       
    71 	/* Initialize tvar to FillValue */
       
    72 	for (i = 0; i < (d1 * d2); i++)
       
    73 		tvar[i] = (double)FillValue;
       
    74 
       
    75 
       
    76 	start[0] = tslice;
       
    77 	count[0] = 1;
       
    78 	start[2] = 0;
       
    79 	count[2] = d1;
       
    80 	start[3] = 0;
       
    81 	count[3] = d2;
       
    82 
       
    83 	/* March through the levels, totalling as we go */
       
    84 	for (j = 0; j < nlevels; j++) {
       
    85 		start[1] = j;
       
    86 		count[1] = 1;
       
    87 		wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
       
    88 		for (i = 0; i < (d1 * d2); i++) {
       
    89 			if (tvar[i] == (double)FillValue)
       
    90 				tvar[i] = (double)var[i];
       
    91 			else
       
    92 				tvar[i] += (double)var[i];
       
    93 		}
       
    94 	}
       
    95 
       
    96 	return;
       
    97 }
       
    98 
       
    99 void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)
       
   100 {
       
   101 	int i;
       
   102 	float **svars;
       
   103 
       
   104 	if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {
       
   105 		perror("alloc_svars:svarsp");
       
   106 		exit(2);
       
   107 	}
       
   108 	svars = *svarsp;
       
   109 	for (i = 0; i < nsvars; i++) {
       
   110 		if (!(svars[i] = (float *)malloc(d1 * d2 * sizeof(float)))) {
       
   111 			perror("alloc_svars:svars[i]");
       
   112 			exit(2);
       
   113 		}
       
   114 	}
       
   115 	return;
       
   116 }
       
   117 
       
   118 void free_svars(size_t nsvars, float **svars)
       
   119 {
       
   120 	int i;
       
   121 
       
   122 	for (i = 0; i < nsvars; i++)
       
   123 		free(svars[i]);
       
   124 	free(svars);
       
   125 	return;
       
   126 }
       
   127 
       
   128 void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,
       
   129 	float *var)
       
   130 {
       
   131 	int varid;
       
   132 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
       
   133 
       
   134 	start[0] = tslice;
       
   135 	count[0] = 1;
       
   136 	start[1] = 0;
       
   137 	count[1] = d1;
       
   138 	start[2] = 0;
       
   139 	count[2] = d2;
       
   140 
       
   141 	/* Get variable id */
       
   142 	wrap_nc(nc_inq_varid(ncid, name, &varid));
       
   143 	/* Assume the correct size and read the variable */
       
   144 	wrap_nc(nc_get_vara_float(ncid, varid, start, count, var));
       
   145 
       
   146 	return;
       
   147 }
       
   148 
       
   149 void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,
       
   150 	int num, double *tvar)
       
   151 {
       
   152 	int i;
       
   153 	size_t nsvars;
       
   154 	float **svars; /* source variables */
       
   155 	double denom;
       
   156 
       
   157 	/* Initialize tvar to FillValue */
       
   158 	for (i = 0; i < (d1 * d2); i++)
       
   159 		tvar[i] = (double)combo_FillValues[num];
       
   160 
       
   161 	switch (num) {
       
   162 		case 0: /* LATENT */
       
   163 			nsvars = 3;
       
   164 			alloc_svars(nsvars, d1, d2, &svars);
       
   165 			/* read timeslices */
       
   166 			read_timeslice(ncid, d1, d2, tslice, "FCTR", svars[0]);
       
   167 			read_timeslice(ncid, d1, d2, tslice, "FCEV", svars[1]);
       
   168 			read_timeslice(ncid, d1, d2, tslice, "FGEV", svars[2]);
       
   169 			/* compute new variable */
       
   170 			for (i = 0; i < (d1 * d2); i++) {
       
   171 				if (svars[0][i] != combo_FillValues[num] &&
       
   172 					svars[1][i] != combo_FillValues[num] &&
       
   173 					svars[2][i] != combo_FillValues[num]) {
       
   174 					tvar[i] = (double)svars[0][i]
       
   175 						+ (double)svars[1][i]
       
   176 						+ (double)svars[2][i];
       
   177 				}
       
   178 			}
       
   179 			free_svars(nsvars, svars);
       
   180 			break;
       
   181 		case 1: /* ET */
       
   182 			nsvars = 3;
       
   183 			alloc_svars(nsvars, d1, d2, &svars);
       
   184 			/* read timeslices */
       
   185 			read_timeslice(ncid, d1, d2, tslice, "QVEGE", svars[0]);
       
   186 			read_timeslice(ncid, d1, d2, tslice, "QVEGT", svars[1]);
       
   187 			read_timeslice(ncid, d1, d2, tslice, "QSOIL", svars[2]);
       
   188 			/* compute new variable */
       
   189 			for (i = 0; i < (d1 * d2); i++) {
       
   190 				if (svars[0][i] != combo_FillValues[num] &&
       
   191 					svars[1][i] != combo_FillValues[num] &&
       
   192 					svars[2][i] != combo_FillValues[num]) {
       
   193 					tvar[i] = (double)svars[0][i]
       
   194 						+ (double)svars[1][i]
       
   195 						+ (double)svars[2][i];
       
   196 				}
       
   197 			}
       
   198 			free_svars(nsvars, svars);
       
   199 			break;
       
   200 		case 2: /* ALL_SKY_ALBEDO */
       
   201 			nsvars = 2;
       
   202 			alloc_svars(nsvars, d1, d2, &svars);
       
   203 			/* read timeslices */
       
   204 			read_timeslice(ncid, d1, d2, tslice, "FSR", svars[0]);
       
   205 			read_timeslice(ncid, d1, d2, tslice, "FSDS", svars[1]);
       
   206 			/* compute new variable */
       
   207 			for (i = 0; i < (d1 * d2); i++) {
       
   208 				if (svars[0][i] != combo_FillValues[num] &&
       
   209 					svars[1][i] != combo_FillValues[num] &&
       
   210 					(svars[1][i] > 1.e-35 || svars[1][i] < -1.e-35)) {
       
   211 					tvar[i] = (double)svars[0][i]
       
   212 						/ (double)svars[1][i];
       
   213 				}
       
   214 			}
       
   215 			free_svars(nsvars, svars);
       
   216 			break;
       
   217 		case 3: /* BLACK_SKY_ALBEDO */
       
   218 			nsvars = 4;
       
   219 			alloc_svars(nsvars, d1, d2, &svars);
       
   220 			/* read timeslices */
       
   221 			read_timeslice(ncid, d1, d2, tslice, "FSRNDLN", svars[0]);
       
   222 			read_timeslice(ncid, d1, d2, tslice, "FSRVDLN", svars[1]);
       
   223 			read_timeslice(ncid, d1, d2, tslice, "FSDSNDLN", svars[2]);
       
   224 			read_timeslice(ncid, d1, d2, tslice, "FSDSVDLN", svars[3]);
       
   225 			/* compute new variable */
       
   226 			for (i = 0; i < (d1 * d2); i++) {
       
   227 				if (svars[0][i] != combo_FillValues[num] &&
       
   228 					svars[1][i] != combo_FillValues[num] &&
       
   229 					svars[2][i] != combo_FillValues[num] &&
       
   230 					svars[3][i] != combo_FillValues[num]) {
       
   231 					denom = (double)svars[2][i] + (double)svars[3][i];
       
   232 					if (denom > 1.e-35 || denom < -1.e-35) {
       
   233 						tvar[i] = ((double)svars[0][i]
       
   234 							+ (double)svars[1][i]) /
       
   235 							((double)svars[2][i]
       
   236 							+ (double)svars[3][i]);
       
   237 					}
       
   238 				}
       
   239 			}
       
   240 			free_svars(nsvars, svars);
       
   241 			break;
       
   242 		case 4: /* NETRAD */
       
   243 			nsvars = 2;
       
   244 			alloc_svars(nsvars, d1, d2, &svars);
       
   245 			/* read timeslices */
       
   246 			read_timeslice(ncid, d1, d2, tslice, "FSA", svars[0]);
       
   247 			read_timeslice(ncid, d1, d2, tslice, "FIRA", svars[1]);
       
   248 			/* compute new variable */
       
   249 			for (i = 0; i < (d1 * d2); i++) {
       
   250 				if (svars[0][i] != combo_FillValues[num] &&
       
   251 					svars[1][i] != combo_FillValues[num]) {
       
   252 					tvar[i] = (double)svars[0][i]
       
   253 						- (double)svars[1][i];
       
   254 				}
       
   255 			}
       
   256 			free_svars(nsvars, svars);
       
   257 			break;
       
   258 		default:
       
   259 			fprintf(stderr, "Unknown combination variable %d!\n", num);
       
   260 			exit(3);
       
   261 	}
       
   262 	
       
   263 	return;
       
   264 }
       
   265 
       
   266 void write_timeslice(int ncid, int varid, size_t d1, size_t d2,
       
   267 	size_t tslice, float *var)
       
   268 {
       
   269 	size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
       
   270 
       
   271 	start[0] = tslice;
       
   272 	count[0] = 1;
       
   273 	start[1] = 0;
       
   274 	count[1] = d1;
       
   275 	start[2] = 0;
       
   276 	count[2] = d2;
       
   277 
       
   278 	wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));
       
   279 
       
   280 	return;
       
   281 }
       
   282 
       
   283 void usage(char *program)
       
   284 {
       
   285 	fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",
       
   286 		program);
       
   287 	fprintf(stderr, "Requires at least one history output file name\n");
       
   288 	exit(5);
       
   289 }
       
   290 
       
   291 int main(int argc, char **argv)
       
   292 {
       
   293 	char **ifnames;
       
   294 	int i, j, k, nifnames;
       
   295 	int ncid, ngdims, nvars, ngatts, unlimdimid, varid, ndims, natts,
       
   296 		dimids[NC_MAX_VAR_DIMS], new_ndims, new_dimids[NC_MAX_VAR_DIMS],
       
   297 		new_varid, lat_dimid, lon_dimid;
       
   298 	nc_type xtype, atype;
       
   299 	char name[NC_MAX_NAME+1], new_name[NC_MAX_NAME+1],
       
   300 		att_name[NC_MAX_NAME+1], *longname, *new_longname;
       
   301 	size_t tlen, ts, alen, nlevels, d1, d2;
       
   302 	float FillValue;
       
   303 	char FillFlag;
       
   304 	float *var;
       
   305 	double *tvar;
       
   306 
       
   307 #ifdef DEBUG
       
   308 	printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);
       
   309 	printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);
       
   310 #endif /* DEBUG */
       
   311 
       
   312 	/* Require at least one parameter (the file on which to work) */
       
   313 	if (argc < 2)
       
   314 		usage(argv[0]);
       
   315 
       
   316 	/* Put the list of filenames in an array of strings */
       
   317 	if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {
       
   318 		perror("ifnames");
       
   319 		exit(2);
       
   320 	}
       
   321 	nifnames = 0;
       
   322 	for (i = 1; i < argc; i++) {
       
   323 		if (!(ifnames[nifnames++] = strdup(argv[i]))) {
       
   324 			perror("ifnames[nifnames++]");
       
   325 			exit(2);
       
   326 		}
       
   327 	}
       
   328 
       
   329 	/*
       
   330 	 * March through all input files adding new variables for totals
       
   331 	 * across levels.
       
   332 	 */
       
   333 	for (i = 0; i < nifnames; i++) {
       
   334 		printf("Processing %s\n", ifnames[i]);
       
   335 		/* Open file */
       
   336 		wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));
       
   337 		/*
       
   338 		 * Inquire about number of dimensions, variables, global
       
   339 		 * attributes and the ID of the unlimited dimension
       
   340 		 */
       
   341 		wrap_nc(nc_inq(ncid, &ngdims, &nvars, &ngatts, &unlimdimid));
       
   342 		wrap_nc(nc_inq_dim(ncid, unlimdimid, name, &tlen));
       
   343 		if (strcmp(name, time_name)) {
       
   344 			fprintf(stderr, "%s is no the unlimited dimension for file %s!\n", time_name, ifnames[i]);
       
   345 			exit(4);
       
   346 		}
       
   347 		/* Retrieve lat and lon dimension IDs */
       
   348 		wrap_nc(nc_inq_dimid(ncid, lat_name, &lat_dimid));
       
   349 		wrap_nc(nc_inq_dimid(ncid, lon_name, &lon_dimid));
       
   350 		/* Put file into define mode */
       
   351 		wrap_nc(nc_redef(ncid));
       
   352 		/*
       
   353 		 * Define all of the new variables first to improve
       
   354 		 * performance.  This means some calls are made twice.
       
   355 		 */
       
   356 		/* First, add the new total variables */
       
   357 		for (j = 0; j < nmtotal_vars; j++) {
       
   358 			/* Figure out source variable information */
       
   359 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
       
   360 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
       
   361 			/* Make sure it has the correct number of dimensions */
       
   362 			if (ndims != 4) {
       
   363 				fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);
       
   364 				exit(4);
       
   365 			}
       
   366 			/* Make sure time is the first dimension */
       
   367 			if (dimids[0] != unlimdimid) {
       
   368 				fprintf(stderr, "First dimension of variable %s is not %s!\n", name, time_name);
       
   369 				exit(4);
       
   370 			}
       
   371 			/* Make sure it is a float type */
       
   372 			if (xtype != NC_FLOAT) {
       
   373 				fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");
       
   374 				exit(4);
       
   375 			}
       
   376 			/* Set dimensions for new variable */
       
   377 			new_ndims = 3;
       
   378 			new_dimids[0] = unlimdimid;
       
   379 			new_dimids[1] = dimids[2];
       
   380 			new_dimids[2] = dimids[3];
       
   381 			/* Define new variable */
       
   382 			sprintf(new_name, "%s%s", varname_total_prefix, name);
       
   383 			printf("\tAdding variable %s\n", new_name);
       
   384 			wrap_nc(nc_def_var(ncid, new_name, xtype, new_ndims,
       
   385 				new_dimids, &new_varid));
       
   386 			/* Copy attributes from source variable to new one */
       
   387 			for (k = 0; k < natts; k++) {
       
   388 				wrap_nc(nc_inq_attname(ncid, varid, k, att_name));
       
   389 				if (!strcmp(att_name, "long_name")) {
       
   390 					wrap_nc(nc_inq_att(ncid, varid, att_name, &atype, &alen));
       
   391 					if (!(longname = (char *)malloc((alen+1)*sizeof(char)))) {
       
   392 						perror("longname");
       
   393 						exit(2);
       
   394 					}
       
   395 					wrap_nc(nc_get_att_text(ncid, varid, att_name, longname));
       
   396 					longname[alen] = '\0';
       
   397 					if (!(new_longname = (char *)malloc((strlen(longname_total_prefix)+alen+1)*sizeof(char)))) {
       
   398 						perror("new_longname");
       
   399 						exit(2);
       
   400 					}
       
   401 					sprintf(new_longname, "%s%s", longname_total_prefix, longname);
       
   402 					wrap_nc(nc_put_att_text(ncid, new_varid, att_name, strlen(new_longname), new_longname));
       
   403 					free(new_longname);
       
   404 					free(longname);
       
   405 				}
       
   406 				else
       
   407 					wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));
       
   408 			}
       
   409 		}
       
   410 		/* Second, add the new total variables */
       
   411 		for (j = 0; j < nmcombo_vars; j++) {
       
   412 			/* Set dimensions for new variable */
       
   413 			new_ndims = 3;
       
   414 			new_dimids[0] = unlimdimid;
       
   415 			new_dimids[1] = lat_dimid;
       
   416 			new_dimids[2] = lon_dimid;
       
   417 			/* Define new variable */
       
   418 			printf("\tAdding variable %s\n", combo_vars[j]);
       
   419 			wrap_nc(nc_def_var(ncid, combo_vars[j], NC_FLOAT,
       
   420 				new_ndims, new_dimids, &new_varid));
       
   421 			/* Add attributes */
       
   422 			wrap_nc(nc_put_att_text(ncid, new_varid, long_name_name,
       
   423 				strlen(combo_long_names[j]),
       
   424 				combo_long_names[j]));
       
   425 			wrap_nc(nc_put_att_text(ncid, new_varid, units_name,
       
   426 				strlen(combo_units[j]), combo_units[j]));
       
   427 			wrap_nc(nc_put_att_text(ncid, new_varid,
       
   428 				cell_method_name, strlen(cell_method),
       
   429 				cell_method));
       
   430 			wrap_nc(nc_put_att_float(ncid, new_varid,
       
   431 				FillValue_name, NC_FLOAT, 1,
       
   432 				&combo_FillValues[j]));
       
   433 			wrap_nc(nc_put_att_float(ncid, new_varid,
       
   434 				missing_value_name, NC_FLOAT, 1,
       
   435 				&combo_FillValues[j]));
       
   436 		}
       
   437 		/* Leave define mode */
       
   438 		wrap_nc(nc_enddef(ncid));
       
   439 		/* Now actually compute and write out the new total variables */
       
   440 		for (j = 0; j < nmtotal_vars; j++) {
       
   441 			/* Figure out source variable information */
       
   442 			wrap_nc(nc_inq_varid(ncid, total_vars[j], &varid));
       
   443 			wrap_nc(nc_inq_var(ncid, varid, name, &xtype, &ndims, dimids, &natts));
       
   444 			/* Check for _FillValue */
       
   445 			FillFlag = 0;
       
   446 			FillValue = 0.;
       
   447 			if (nc_inq_att(ncid, varid, FillValue_name, &atype, &alen) == NC_NOERR) {
       
   448 				if (atype == NC_FLOAT && alen) {
       
   449 					wrap_nc(nc_get_att_float(ncid, varid,
       
   450 						FillValue_name, &FillValue));
       
   451 					FillFlag = 1;
       
   452 				}
       
   453 			}
       
   454 			/* Set dimensions for new variable */
       
   455 			new_ndims = 3;
       
   456 			new_dimids[0] = unlimdimid;
       
   457 			new_dimids[1] = dimids[2];
       
   458 			new_dimids[2] = dimids[3];
       
   459 
       
   460 			sprintf(new_name, "%s%s", varname_total_prefix, name);
       
   461 			printf("\tComputing and writing %s\n", new_name);
       
   462 			/* Retrieve the new varid again */
       
   463 			wrap_nc(nc_inq_varid(ncid, new_name, &new_varid));
       
   464 			/*
       
   465 			 * Figure out dimensions and total size
       
   466 			 * for a single time slice
       
   467 			 */
       
   468 			wrap_nc(nc_inq_dimlen(ncid, dimids[1], &nlevels));
       
   469 			wrap_nc(nc_inq_dimlen(ncid, dimids[2], &d1));
       
   470 			wrap_nc(nc_inq_dimlen(ncid, dimids[3], &d2));
       
   471 
       
   472 			/*
       
   473 			 * Allocate space for reading in time slice and for
       
   474 			 * the sum
       
   475 			 */
       
   476 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
       
   477 				perror("var");
       
   478 				exit(2);
       
   479 			}
       
   480 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
       
   481 				perror("tvar");
       
   482 				exit(2);
       
   483 			}
       
   484 
       
   485 			/* Read timeslices, total, and write out new timeslices */
       
   486 			for (ts = 0; ts < tlen; ts++) {
       
   487 				read_total_timeslice(ncid, varid,
       
   488 					nlevels, d1, d2, FillFlag, FillValue,
       
   489 					ts, var, tvar);
       
   490 				for (k = 0; k < (d1*d2); k++)
       
   491 					var[k] = (float)tvar[k];
       
   492 				write_timeslice(ncid, new_varid, d1, d2,
       
   493 					ts, var);
       
   494 			}
       
   495 
       
   496 			free(var);
       
   497 			free(tvar);
       
   498 		}
       
   499 		/* Now actually compute and write out the new combo variables */
       
   500 		for (j = 0; j < nmcombo_vars; j++) {
       
   501 			/* Set dimensions for new variable */
       
   502 			new_ndims = 3;
       
   503 			new_dimids[0] = unlimdimid;
       
   504 			new_dimids[1] = lat_dimid;
       
   505 			new_dimids[2] = lon_dimid;
       
   506 
       
   507 			printf("\tComputing and writing %s\n", combo_vars[j]);
       
   508 			/* Retrieve the new varid again */
       
   509 			wrap_nc(nc_inq_varid(ncid, combo_vars[j], &new_varid));
       
   510 			/*
       
   511 			 * Figure out dimensions and total size
       
   512 			 * for a single time slice
       
   513 			 */
       
   514 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));
       
   515 			wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));
       
   516 
       
   517 			/*
       
   518 			 * Allocate space for reading in time slice and for
       
   519 			 * the sum
       
   520 			 */
       
   521 			if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
       
   522 				perror("var");
       
   523 				exit(2);
       
   524 			}
       
   525 			if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
       
   526 				perror("tvar");
       
   527 				exit(2);
       
   528 			}
       
   529 
       
   530 			/* Read timeslices, compute new variable, and write */
       
   531 			for (ts = 0; ts < tlen; ts++) {
       
   532 				read_compute_timeslice(ncid, d1, d2, ts,
       
   533 					j, tvar);
       
   534 				for (k = 0; k < (d1*d2); k++)
       
   535 					var[k] = (float)tvar[k];
       
   536 				write_timeslice(ncid, new_varid, d1, d2,
       
   537 					ts, var);
       
   538 			}
       
   539 
       
   540 			free(var);
       
   541 			free(tvar);
       
   542 		}
       
   543 		/* Close file */
       
   544 		wrap_nc(nc_close(ncid));
       
   545 	}
       
   546 
       
   547 	return 0;
       
   548 }