Fixed h1_summary and h1_summary2 to correctly construct time_bounds values.
h1_summary and h1_summary2 previously used the time_bounds from the last time
value from the last input file when summarizing, suggesting that the field
values were appropriate only over that short time range instead of the complete
time period over which the statistics were calculated. C-LAMP Experiment 1
runs used the previous code, so the time_bounds were incorrect in the
statistical summaries produced. C-LAMP Experiment 2 runs will use this new
code for production of statistical summaries.
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[]
13 static char *total_vars[] = {
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 */
25 static char *combo_long_names[] = {
28 "surface all-sky albedo",
29 "surface black-sky albedo",
32 static char *combo_units[] = {
39 static float combo_FillValues[] = {
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";
54 void wrap_nc(int status)
56 if (status != NC_NOERR) {
57 fprintf(stderr, "netCDF error: %s\n", nc_strerror(status));
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,
69 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
71 /* Initialize tvar to FillValue */
72 for (i = 0; i < (d1 * d2); i++)
73 tvar[i] = (double)FillValue;
83 /* March through the levels, totalling as we go */
84 for (j = 0; j < nlevels; j++) {
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];
92 tvar[i] += (double)var[i];
99 void alloc_svars(size_t nsvars, size_t d1, size_t d2, float ***svarsp)
104 if (!(*svarsp = (float **)malloc(nsvars * sizeof(float **)))) {
105 perror("alloc_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]");
118 void free_svars(size_t nsvars, float **svars)
122 for (i = 0; i < nsvars; i++)
128 void read_timeslice(int ncid, size_t d1, size_t d2, size_t tslice, char *name,
132 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
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));
149 void read_compute_timeslice(int ncid, size_t d1, size_t d2, size_t tslice,
150 int num, double *tvar)
154 float **svars; /* source variables */
157 /* Initialize tvar to FillValue */
158 for (i = 0; i < (d1 * d2); i++)
159 tvar[i] = (double)combo_FillValues[num];
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];
179 free_svars(nsvars, svars);
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];
198 free_svars(nsvars, svars);
200 case 2: /* ALL_SKY_ALBEDO */
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];
215 free_svars(nsvars, svars);
217 case 3: /* BLACK_SKY_ALBEDO */
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]) /
236 + (double)svars[3][i]);
240 free_svars(nsvars, svars);
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];
256 free_svars(nsvars, svars);
259 fprintf(stderr, "Unknown combination variable %d!\n", num);
266 void write_timeslice(int ncid, int varid, size_t d1, size_t d2,
267 size_t tslice, float *var)
269 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS];
278 wrap_nc(nc_put_vara_float(ncid, varid, start, count, var));
283 void usage(char *program)
285 fprintf(stderr, "Usage: %s hist_file1.nc [ hist_file2.nc ... ]\n",
287 fprintf(stderr, "Requires at least one history output file name\n");
291 int main(int argc, char **argv)
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;
308 printf("Number of total variables (nmtotal_vars) = %d\n", nmtotal_vars);
309 printf("Number of combination variables (nmcombo_vars) = %d\n", nmcombo_vars);
312 /* Require at least one parameter (the file on which to work) */
316 /* Put the list of filenames in an array of strings */
317 if (!(ifnames = (char **)malloc((argc - 1) * sizeof(char *)))) {
322 for (i = 1; i < argc; i++) {
323 if (!(ifnames[nifnames++] = strdup(argv[i]))) {
324 perror("ifnames[nifnames++]");
330 * March through all input files adding new variables for totals
333 for (i = 0; i < nifnames; i++) {
334 printf("Processing %s\n", ifnames[i]);
336 wrap_nc(nc_open(ifnames[i], NC_WRITE, &ncid));
338 * Inquire about number of dimensions, variables, global
339 * attributes and the ID of the unlimited dimension
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]);
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));
353 * Define all of the new variables first to improve
354 * performance. This means some calls are made twice.
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 */
363 fprintf(stderr, "Variable %s has %d dimensions!\n", name, ndims);
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);
371 /* Make sure it is a float type */
372 if (xtype != NC_FLOAT) {
373 fprintf(stderr, "Only NC_FLOAT type is supported presently!\n");
376 /* Set dimensions for new variable */
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)))) {
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");
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));
407 wrap_nc(nc_copy_att(ncid, varid, att_name, ncid, new_varid));
410 /* Second, add the new total variables */
411 for (j = 0; j < nmcombo_vars; j++) {
412 /* Set dimensions for new variable */
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));
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),
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]));
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 */
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));
454 /* Set dimensions for new variable */
456 new_dimids[0] = unlimdimid;
457 new_dimids[1] = dimids[2];
458 new_dimids[2] = dimids[3];
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));
465 * Figure out dimensions and total size
466 * for a single time slice
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));
473 * Allocate space for reading in time slice and for
476 if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
480 if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
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,
490 for (k = 0; k < (d1*d2); k++)
491 var[k] = (float)tvar[k];
492 write_timeslice(ncid, new_varid, d1, d2,
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 */
503 new_dimids[0] = unlimdimid;
504 new_dimids[1] = lat_dimid;
505 new_dimids[2] = lon_dimid;
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));
511 * Figure out dimensions and total size
512 * for a single time slice
514 wrap_nc(nc_inq_dimlen(ncid, new_dimids[1], &d1));
515 wrap_nc(nc_inq_dimlen(ncid, new_dimids[2], &d2));
518 * Allocate space for reading in time slice and for
521 if (!(var = (float *)malloc((d1*d2) * sizeof(float)))) {
525 if (!(tvar = (double *)malloc((d1*d2) * sizeof(double)))) {
530 /* Read timeslices, compute new variable, and write */
531 for (ts = 0; ts < tlen; ts++) {
532 read_compute_timeslice(ncid, d1, d2, ts,
534 for (k = 0; k < (d1*d2); k++)
535 var[k] = (float)tvar[k];
536 write_timeslice(ncid, new_varid, d1, d2,
544 wrap_nc(nc_close(ncid));