|
1 #include <stdio.h> |
|
2 #include <unistd.h> |
|
3 #include <stdlib.h> |
|
4 #include <string.h> |
|
5 #include <math.h> |
|
6 #include "netcdf.h" |
|
7 |
|
8 /* |
|
9 * Loop through one month of hourly history tapes from the Community Land Model |
|
10 * and generate monthly statistics (means and standard deviations) of fields |
|
11 * by hour of day. |
|
12 */ |
|
13 |
|
14 #define MEAN_TYPE 0 |
|
15 #define STDDEV_TYPE 1 |
|
16 |
|
17 #define HOURS_PER_DAY 24 |
|
18 #define MIN_PER_HOUR 60 |
|
19 #define SEC_PER_MIN 60 |
|
20 #define SEC_PER_HOUR (MIN_PER_HOUR * SEC_PER_MIN) |
|
21 |
|
22 static char *metadata_vars[] = { |
|
23 "lon", |
|
24 "lat", |
|
25 "lonrof", |
|
26 "latrof", |
|
27 "time", |
|
28 "hour", /* new metadata variable to be added to output files */ |
|
29 "levsoi", |
|
30 "levlak", |
|
31 "edgen", |
|
32 "edgee", |
|
33 "edges", |
|
34 "edgew", |
|
35 "longxy", |
|
36 "latixy", |
|
37 "area", |
|
38 "areaupsc", |
|
39 "topo", |
|
40 "topodnsc", |
|
41 "landfrac", |
|
42 "landmask", |
|
43 "pftmask", |
|
44 "indxupsc", |
|
45 "mcdate", |
|
46 "mcsec", |
|
47 "mdcur", |
|
48 "mscur", |
|
49 "nstep", |
|
50 "time_bounds", |
|
51 "date_written", |
|
52 "time_written" |
|
53 }; |
|
54 |
|
55 struct text_att { |
|
56 char *name; |
|
57 char *value; |
|
58 struct text_att *next; |
|
59 struct text_att *prev; |
|
60 }; |
|
61 |
|
62 struct dim { |
|
63 int dimid; |
|
64 char *name; |
|
65 size_t len; |
|
66 struct dim *next; |
|
67 struct dim *prev; |
|
68 }; |
|
69 |
|
70 struct var { |
|
71 int ncvarid; |
|
72 char *name; |
|
73 nc_type nctype; |
|
74 int ndims; |
|
75 int *dimids; |
|
76 int natts; |
|
77 char metadata; |
|
78 char FillFlag; |
|
79 float FillValue; |
|
80 struct var *next; |
|
81 struct var *prev; |
|
82 }; |
|
83 |
|
84 static char *time_name = "time"; |
|
85 static char *mcsec_name = "mcsec"; |
|
86 static char *history_name = "history"; |
|
87 static char *nsamples_name = "nsamples"; |
|
88 static char *hour_name = "hour", *hour_long_name = "hour of day", |
|
89 *hour_units = "hour"; |
|
90 static char *cell_method_name = "cell_method"; |
|
91 /* input stuff */ |
|
92 static int nmvars = sizeof(metadata_vars)/sizeof(*metadata_vars); |
|
93 static int input_ncid, input_ndims, input_nvars, input_ngatts, input_unlimdimid; |
|
94 static struct text_att *input_gatt_head = NULL; |
|
95 static struct dim *input_dim_head = NULL, **input_dim_idx; |
|
96 static struct var *input_var_head = NULL; |
|
97 /* translation stuff */ |
|
98 static int *idid2mdid, *idid2sdid; /* dimension IDs */ |
|
99 /* output stuff */ |
|
100 static int mean_ncid, mean_ndims, mean_nvars, mean_ngatts, mean_unlimdimid; |
|
101 static int mean_hour_dimid; /* special notes */ |
|
102 static struct dim *mean_dim_head = NULL; |
|
103 static struct var *mean_var_head = NULL; |
|
104 static int stddev_ncid, stddev_ndims, stddev_nvars, stddev_ngatts, stddev_unlimdimid; |
|
105 static int stddev_hour_dimid; /* special notes */ |
|
106 static struct dim *stddev_dim_head = NULL; |
|
107 static struct var *stddev_var_head = NULL; |
|
108 |
|
109 char is_metadata(char *name) |
|
110 { |
|
111 int i; |
|
112 |
|
113 for (i = 0; i < nmvars && strcmp(name, metadata_vars[i]); i++); |
|
114 |
|
115 if (i < nmvars) |
|
116 return 1; |
|
117 else |
|
118 return 0; |
|
119 } |
|
120 |
|
121 #if 0 |
|
122 struct dim *dimlist_find_by_name(struct dim *head, char *name) |
|
123 { |
|
124 int i; |
|
125 struct dim *p = NULL; |
|
126 |
|
127 if (head) { |
|
128 node = head; |
|
129 for (i = 0 ; i == 0 || node != head; i++) { |
|
130 if (!strcmp(node->name, name)) |
|
131 p = node; |
|
132 node = node->next; |
|
133 } |
|
134 return p; |
|
135 } |
|
136 else |
|
137 return NULL; |
|
138 } |
|
139 #endif |
|
140 |
|
141 struct var *varlist_find_by_name(struct var *head, char *name) |
|
142 { |
|
143 int i; |
|
144 struct var *node; |
|
145 |
|
146 if (head) { |
|
147 node = head; |
|
148 for (i = 0 ; (i == 0 || node != head) && strcmp(node->name, name); i++) |
|
149 node = node->next; |
|
150 if (i && node == head) |
|
151 return NULL; |
|
152 else |
|
153 return node; |
|
154 } |
|
155 else |
|
156 return NULL; |
|
157 } |
|
158 |
|
159 void gattlist_add_node(struct text_att **headp, char *name, char *value) |
|
160 { |
|
161 struct text_att *head, *node; |
|
162 |
|
163 head = *headp; |
|
164 |
|
165 if (!(node = (struct text_att *)malloc(sizeof(struct text_att)))) { |
|
166 perror("gattlist_add_node"); |
|
167 exit(2); |
|
168 } |
|
169 |
|
170 if (!(node->name = strdup(name))) { |
|
171 perror("gattlist_add_node"); |
|
172 exit(2); |
|
173 } |
|
174 if (!(node->value = strdup(value))) { |
|
175 perror("gattlist_add_node"); |
|
176 exit(2); |
|
177 } |
|
178 |
|
179 if (head == NULL) { |
|
180 /* first one */ |
|
181 *headp = node; |
|
182 node->next = node; |
|
183 node->prev = node; |
|
184 } |
|
185 else { |
|
186 node->next = head; |
|
187 node->prev = head->prev; |
|
188 /* set this after setting node->prev from here */ |
|
189 head->prev = node; |
|
190 /* set this after having set node->prev */ |
|
191 node->prev->next = node; |
|
192 } |
|
193 |
|
194 return; |
|
195 } |
|
196 |
|
197 struct dim *dimlist_add_node(struct dim **headp, int dimid, char *name, size_t len) |
|
198 { |
|
199 struct dim *head, *node; |
|
200 |
|
201 head = *headp; |
|
202 |
|
203 if (!(node = (struct dim *)malloc(sizeof(struct dim)))) { |
|
204 perror("dimlist_add_node"); |
|
205 exit(2); |
|
206 } |
|
207 |
|
208 node->dimid = dimid; |
|
209 if (!(node->name = strdup(name))) { |
|
210 perror("dimlist_add_node"); |
|
211 exit(2); |
|
212 } |
|
213 node->len = len; |
|
214 |
|
215 if (head == NULL) { |
|
216 /* first one */ |
|
217 *headp = node; |
|
218 node->next = node; |
|
219 node->prev = node; |
|
220 } |
|
221 else { |
|
222 node->next = head; |
|
223 node->prev = head->prev; |
|
224 /* set this after setting node->prev from here */ |
|
225 head->prev = node; |
|
226 /* set this after having set node->prev */ |
|
227 node->prev->next = node; |
|
228 } |
|
229 |
|
230 return node; |
|
231 } |
|
232 |
|
233 void varlist_add_node(struct var **headp, int ncvarid, char *name, |
|
234 nc_type nctype, int ndims, int *dimids, int natts, char FillFlag, |
|
235 float FillValue) |
|
236 { |
|
237 int i; |
|
238 struct var *head, *node; |
|
239 |
|
240 head = *headp; |
|
241 |
|
242 if (!(node = (struct var *)malloc(sizeof(struct var)))) { |
|
243 perror("varlist_add_node"); |
|
244 exit(3); |
|
245 } |
|
246 |
|
247 node->ncvarid = ncvarid; |
|
248 if (!(node->name = strdup(name))) { |
|
249 perror("varlist_add_node"); |
|
250 exit(3); |
|
251 } |
|
252 node->nctype = nctype; |
|
253 node->ndims = ndims; |
|
254 if (!(node->dimids = (int *)malloc(ndims * sizeof(int)))) { |
|
255 perror("varlist_add_node"); |
|
256 exit(3); |
|
257 } |
|
258 for (i = 0; i < ndims; i++) node->dimids[i] = dimids[i]; |
|
259 node->natts = natts; |
|
260 node->metadata = is_metadata(name); |
|
261 node->FillFlag = FillFlag; |
|
262 node->FillValue = FillValue; |
|
263 |
|
264 if (head == NULL) { |
|
265 /* first one */ |
|
266 *headp = node; |
|
267 node->next = node; |
|
268 node->prev = node; |
|
269 } |
|
270 else { |
|
271 node->next = head; |
|
272 node->prev = head->prev; |
|
273 /* set this after setting node->prev from here */ |
|
274 head->prev = node; |
|
275 /* set this after having set node->prev */ |
|
276 node->prev->next = node; |
|
277 } |
|
278 |
|
279 return; |
|
280 } |
|
281 |
|
282 void varlist_print(struct var *headp) |
|
283 { |
|
284 int i, j; |
|
285 struct var *node; |
|
286 |
|
287 printf("Beginning of Variable List\n"); |
|
288 |
|
289 if (headp) { |
|
290 node = headp; |
|
291 for (j = 0; j == 0 || node != headp; j++) { |
|
292 printf("Variable %d (ptr=%ld):\n", j, (long)node); |
|
293 printf("\tncvarid = %d\n", node->ncvarid); |
|
294 printf("\tname = %s\n", node->name); |
|
295 printf("\tnctype = %d\n", node->nctype); |
|
296 printf("\tndims = %d\n", node->ndims); |
|
297 printf("\tdimids ="); |
|
298 for (i = 0; i < node->ndims; i++) |
|
299 printf(" %d", node->dimids[i]); |
|
300 printf("\n"); |
|
301 printf("\tnatts = %d\n", node->natts); |
|
302 printf("\tmetadata = %d\n", (int)node->metadata); |
|
303 printf("\tnext = %ld\n", (long)node->next); |
|
304 printf("\tprev = %ld\n", (long)node->prev); |
|
305 node = node->next; |
|
306 } |
|
307 } |
|
308 else { |
|
309 printf("\tList undefined\n"); |
|
310 } |
|
311 |
|
312 printf("End of Variable List\n"); |
|
313 |
|
314 return; |
|
315 } |
|
316 |
|
317 void wrap_nc(int status) |
|
318 { |
|
319 if (status != NC_NOERR) { |
|
320 fprintf(stderr, "netCDF error: %s\n", nc_strerror(status)); |
|
321 exit(1); |
|
322 } |
|
323 |
|
324 return; |
|
325 } |
|
326 |
|
327 void get_dimensions(int ncid, int ndims, struct dim **dim_headp, struct dim ***in_dim_idxp) |
|
328 { |
|
329 int i; |
|
330 char name[NC_MAX_NAME+1]; |
|
331 size_t len; |
|
332 struct dim **in_dim_idx; |
|
333 |
|
334 if (!(*in_dim_idxp = (struct dim **)malloc(ndims * sizeof(struct dim *)))) { |
|
335 perror("*in_dim_idxp"); |
|
336 exit(3); |
|
337 } |
|
338 in_dim_idx = *in_dim_idxp; |
|
339 |
|
340 for (i = 0; i < ndims; i++) { |
|
341 wrap_nc(nc_inq_dim(ncid, i, name, &len)); |
|
342 in_dim_idx[i] = dimlist_add_node(dim_headp, i, name, len); |
|
343 } |
|
344 |
|
345 return; |
|
346 } |
|
347 void get_gatts(int ncid, int ngatts, struct text_att **gatt_headp) |
|
348 { |
|
349 int i; |
|
350 char name[NC_MAX_NAME+1], *value; |
|
351 nc_type xtype; |
|
352 size_t len; |
|
353 |
|
354 for (i = 0; i < ngatts; i++) { |
|
355 wrap_nc(nc_inq_attname(ncid, NC_GLOBAL, i, name)); |
|
356 wrap_nc(nc_inq_att(ncid, NC_GLOBAL, name, &xtype, &len)); |
|
357 if (xtype != NC_CHAR) { |
|
358 fprintf(stderr, "Global attribute %s is not of type NC_CHAR\n", name); |
|
359 exit(2); |
|
360 } |
|
361 if (!(value = (char *)malloc((len+1)*sizeof(char)))) { |
|
362 perror("get_gatts: value"); |
|
363 exit(3); |
|
364 } |
|
365 wrap_nc(nc_get_att_text(ncid, NC_GLOBAL, name, value)); |
|
366 value[(len+1-1)] = '\0'; |
|
367 gattlist_add_node(gatt_headp, name, value); |
|
368 free(value); /* because gattlist_add_node() duplicates it */ |
|
369 } |
|
370 |
|
371 return; |
|
372 } |
|
373 |
|
374 void get_vars(int ncid, int nvars, struct var **var_headp) |
|
375 { |
|
376 int i, ndims, dimids[NC_MAX_VAR_DIMS], natts; |
|
377 size_t len; |
|
378 float FillValue; |
|
379 char name[NC_MAX_NAME+1], *fill_att_name = "_FillValue", FillFlag; |
|
380 nc_type xtype, atype; |
|
381 |
|
382 for (i = 0; i < nvars; i++) { |
|
383 wrap_nc(nc_inq_var(ncid, i, name, &xtype, &ndims, dimids, |
|
384 &natts)); |
|
385 /* Check for _FillValue */ |
|
386 FillFlag = 0; |
|
387 FillValue = 0.; |
|
388 if (nc_inq_att(ncid, i, fill_att_name, &atype, &len) |
|
389 == NC_NOERR) { |
|
390 if (atype == NC_FLOAT && len) { |
|
391 wrap_nc(nc_get_att_float(ncid, i, |
|
392 fill_att_name, &FillValue)); |
|
393 FillFlag = 1; |
|
394 } |
|
395 } |
|
396 varlist_add_node(var_headp, i, name, xtype, ndims, dimids, |
|
397 natts, FillFlag, FillValue); |
|
398 } |
|
399 |
|
400 return; |
|
401 } |
|
402 |
|
403 int put_dimensions(struct dim *in_dim_head, int in_ndims, int in_unlimdimid, |
|
404 size_t nsamples, int **in2outp, int out_ncid, |
|
405 struct dim **out_dim_headp, int *out_unlimdimidp, int *out_hour_dimidp) |
|
406 { |
|
407 /* |
|
408 * Define dimensions on new files and build translation tables between |
|
409 * dimension IDs. |
|
410 */ |
|
411 int j, dimid, ndims, *in2out; |
|
412 size_t len; |
|
413 struct dim *dnode; |
|
414 |
|
415 ndims = 0; |
|
416 |
|
417 if (!(*in2outp = (int *)malloc((in_ndims+1)*sizeof(int)))) { |
|
418 perror("put_dimensions: in2outp"); |
|
419 exit(3); |
|
420 } |
|
421 in2out = *in2outp; |
|
422 |
|
423 if (in_dim_head) { |
|
424 dnode = in_dim_head; |
|
425 for (j = 0; j == 0 || dnode != in_dim_head; j++) { |
|
426 if (dnode->dimid == in_unlimdimid) |
|
427 len = NC_UNLIMITED; |
|
428 else |
|
429 len = dnode->len; |
|
430 wrap_nc(nc_def_dim(out_ncid, dnode->name, len, &dimid)); |
|
431 dimlist_add_node(out_dim_headp, dimid, dnode->name, len); |
|
432 ++ndims; |
|
433 in2out[dnode->dimid] = dimid; |
|
434 if (dnode->dimid == in_unlimdimid) |
|
435 *out_unlimdimidp = dimid; |
|
436 /* |
|
437 * Just after the "time" dimension, add the new |
|
438 * "nsamples" and "hour" dimensions. |
|
439 */ |
|
440 if (!strcmp(dnode->name, time_name)) { |
|
441 wrap_nc(nc_def_dim(out_ncid, nsamples_name, nsamples, &dimid)); |
|
442 dimlist_add_node(out_dim_headp, dimid, nsamples_name, nsamples); |
|
443 ++ndims; |
|
444 |
|
445 wrap_nc(nc_def_dim(out_ncid, hour_name, HOURS_PER_DAY, &dimid)); |
|
446 dimlist_add_node(out_dim_headp, dimid, hour_name, HOURS_PER_DAY); |
|
447 ++ndims; |
|
448 /* Track hour dimid for out files */ |
|
449 *out_hour_dimidp = dimid; |
|
450 } |
|
451 |
|
452 dnode = dnode->next; |
|
453 } |
|
454 } |
|
455 else { |
|
456 fprintf(stderr, "WARNING: No dimensions defined!\n"); |
|
457 } |
|
458 |
|
459 return ndims; |
|
460 } |
|
461 |
|
462 int put_gatts(struct text_att *in_gatt_head, int out_ncid, char *out_history) |
|
463 { |
|
464 /* |
|
465 * Write out global attributes matching those of the input file. |
|
466 * Change history attribute to the string passed in as an argument. |
|
467 */ |
|
468 int j, hflag = 0, ngatts; |
|
469 char *value; |
|
470 struct text_att *anode; |
|
471 |
|
472 ngatts = 0; |
|
473 |
|
474 if (in_gatt_head) { |
|
475 anode = in_gatt_head; |
|
476 for (j = 0; j == 0 || anode != in_gatt_head; j++) { |
|
477 if (!strcmp(anode->name, history_name)) { |
|
478 value = out_history; |
|
479 ++hflag; |
|
480 } |
|
481 else |
|
482 value = anode->value; |
|
483 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, anode->name, strlen(value), value)); |
|
484 ++ngatts; |
|
485 anode = anode->next; |
|
486 } |
|
487 /* If no history attribute on input, add one to the output */ |
|
488 if (!hflag) { |
|
489 wrap_nc(nc_put_att_text(out_ncid, NC_GLOBAL, history_name, strlen(out_history), out_history)); |
|
490 ++ngatts; |
|
491 } |
|
492 } |
|
493 else { |
|
494 fprintf(stderr, "WARNING: No global attributes defined!\n"); |
|
495 } |
|
496 |
|
497 return ngatts; |
|
498 } |
|
499 |
|
500 int translate_dimids(struct dim **in_dim_idx, char metadata, int in_ndims, int *in_dimids, int *in2out, int *out_dimids, int hour_dimid) |
|
501 { |
|
502 /* |
|
503 * Translate between input dimension IDs and those for the output file. |
|
504 * For normal time-based variables, add a new dimension for hour of |
|
505 * day. For metadata variables, do not add new dimensions, even if |
|
506 * it is time-based. Return (possibly new) number of dimensions. |
|
507 */ |
|
508 int i, ndims; |
|
509 |
|
510 for (i = ndims = 0; i < in_ndims; i++, ndims++) { |
|
511 out_dimids[ndims] = in2out[in_dimids[i]]; |
|
512 if (!strcmp((in_dim_idx[in_dimids[i]])->name, time_name) |
|
513 && !metadata) { |
|
514 ndims++; |
|
515 out_dimids[ndims] = hour_dimid; |
|
516 } |
|
517 } |
|
518 |
|
519 return ndims; |
|
520 } |
|
521 |
|
522 int copy_att(char metadata, int stat_type, int in_natts, int in_ncid, |
|
523 int in_varid, int out_ncid, int out_varid) |
|
524 { |
|
525 /* |
|
526 * Copy the attributes of the input variable to those of the output |
|
527 * variable. If the variable is not a metadata variable, ensure that |
|
528 * the cell_method attribute is set either to "time: mean" or |
|
529 * "time: standard_deviation", depending on the output file type. |
|
530 */ |
|
531 |
|
532 int i, natts; |
|
533 char name[NC_MAX_NAME+1], cmflag = 0; |
|
534 char *cell_methods[] = { "time: mean", "time: standard_deviation" }; |
|
535 |
|
536 for (i = natts = 0; i < in_natts; i++, natts++) { |
|
537 wrap_nc(nc_inq_attname(in_ncid, in_varid, i, name)); |
|
538 if (!strcmp(name, cell_method_name) && !metadata) { |
|
539 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); |
|
540 cmflag = 1; |
|
541 } |
|
542 else |
|
543 wrap_nc(nc_copy_att(in_ncid, in_varid, name, out_ncid, out_varid)); |
|
544 } |
|
545 /* |
|
546 * If no cell_method attribute was specified for a non-metadata |
|
547 * variable on the input file, add the appropriate cell_method anyway |
|
548 * on the output file. |
|
549 */ |
|
550 if (!cmflag && !metadata) { |
|
551 wrap_nc(nc_put_att_text(out_ncid, out_varid, cell_method_name, strlen(cell_methods[stat_type]), cell_methods[stat_type])); |
|
552 natts++; |
|
553 } |
|
554 |
|
555 return natts; |
|
556 } |
|
557 |
|
558 int define_vars(int in_ncid, struct dim **in_dim_idx, struct var *in_var_head, |
|
559 int *in2out, int out_ncid, int hour_dimid, struct var **out_var_headp, |
|
560 int stat_type) |
|
561 { |
|
562 /* |
|
563 * Define variables on output file and copy attributes from input file. |
|
564 * Return number of variables defined. |
|
565 */ |
|
566 int j, varid, nvars, ndims, dimids[NC_MAX_VAR_DIMS], natts; |
|
567 struct var *vnode; |
|
568 |
|
569 nvars = 0; |
|
570 |
|
571 if (in_var_head) { |
|
572 vnode = in_var_head; |
|
573 /* |
|
574 * March through input variables creating (mostly) the same |
|
575 * variables on the output file. |
|
576 */ |
|
577 for (j = 0; j == 0 || vnode != in_var_head; j++) { |
|
578 ndims = translate_dimids(in_dim_idx, vnode->metadata, vnode->ndims, vnode->dimids, in2out, dimids, hour_dimid); |
|
579 wrap_nc(nc_def_var(out_ncid, vnode->name, vnode->nctype, ndims, dimids, &varid)); |
|
580 natts = copy_att(vnode->metadata, stat_type, vnode->natts, in_ncid, vnode->ncvarid, out_ncid, varid); |
|
581 varlist_add_node(out_var_headp, varid, vnode->name, vnode->nctype, ndims, dimids, natts, vnode->FillFlag, vnode->FillValue); |
|
582 ++nvars; |
|
583 /* |
|
584 * Just after the "time" variable, add the new "hour" |
|
585 * variable for hour of day. |
|
586 */ |
|
587 if (!strcmp(vnode->name, time_name)) { |
|
588 ndims = 1; |
|
589 dimids[0] = hour_dimid; |
|
590 wrap_nc(nc_def_var(out_ncid, hour_name, NC_FLOAT, ndims, dimids, &varid)); |
|
591 wrap_nc(nc_put_att_text(out_ncid, varid, "long_name", strlen(hour_long_name), hour_long_name)); |
|
592 wrap_nc(nc_put_att_text(out_ncid, varid, "units", strlen(hour_units), hour_units)); |
|
593 varlist_add_node(out_var_headp, varid, hour_name, NC_FLOAT, ndims, dimids, 2, 0, 0.0); |
|
594 ++nvars; |
|
595 } |
|
596 |
|
597 vnode = vnode->next; |
|
598 } |
|
599 } |
|
600 else { |
|
601 fprintf(stderr, "WARNING: No variables defined!\n"); |
|
602 } |
|
603 |
|
604 return nvars; |
|
605 } |
|
606 |
|
607 void *alloc_var(nc_type nctype, size_t len) |
|
608 { |
|
609 void *val; |
|
610 |
|
611 switch (nctype) { |
|
612 case NC_FLOAT: |
|
613 if (!(val = (float *)malloc((len) * sizeof(float)))) { |
|
614 perror("alloc_var: val"); |
|
615 exit(3); |
|
616 } |
|
617 break; |
|
618 case NC_INT: |
|
619 if (!(val = (int *)malloc((len) * sizeof(int)))) { |
|
620 perror("alloc_var: val"); |
|
621 exit(3); |
|
622 } |
|
623 break; |
|
624 case NC_DOUBLE: |
|
625 if (!(val = (double *)malloc((len) * sizeof(double)))) { |
|
626 perror("alloc_var: val"); |
|
627 exit(3); |
|
628 } |
|
629 break; |
|
630 case NC_CHAR: |
|
631 if (!(val = (char *)malloc(((len)+1) * sizeof(char)))) { |
|
632 perror("alloc_var: val"); |
|
633 exit(3); |
|
634 } |
|
635 break; |
|
636 default: |
|
637 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
638 exit(3); |
|
639 } |
|
640 |
|
641 return val; |
|
642 } |
|
643 |
|
644 void *read_var(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx) |
|
645 { |
|
646 int i; |
|
647 size_t len = (size_t)1; |
|
648 void *val; |
|
649 |
|
650 /* Figure out the total size */ |
|
651 for (i = 0; i < ndims; i++) len *= (dim_idx[dimids[i]])->len; |
|
652 |
|
653 val = alloc_var(nctype, len); |
|
654 |
|
655 switch (nctype) { |
|
656 case NC_FLOAT: |
|
657 wrap_nc(nc_get_var_float(ncid, varid, val)); |
|
658 break; |
|
659 case NC_INT: |
|
660 wrap_nc(nc_get_var_int(ncid, varid, val)); |
|
661 break; |
|
662 case NC_DOUBLE: |
|
663 wrap_nc(nc_get_var_double(ncid, varid, val)); |
|
664 break; |
|
665 case NC_CHAR: |
|
666 wrap_nc(nc_get_var_text(ncid, varid, val)); |
|
667 break; |
|
668 default: |
|
669 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
670 exit(3); |
|
671 } |
|
672 |
|
673 return val; |
|
674 } |
|
675 |
|
676 void *read_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, size_t tslice) |
|
677 { |
|
678 int i; |
|
679 size_t len = (size_t)1, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; |
|
680 void *val; |
|
681 |
|
682 /* Make sure time is really the first dimension */ |
|
683 if (strcmp((dim_idx[dimids[0]])->name, time_name)) { |
|
684 fprintf(stderr, "read_timeslice: %s is not first dimension of variable!\n", time_name); |
|
685 exit(3); |
|
686 } |
|
687 /* |
|
688 * Figure out the total size for the slice (start at second dimension) |
|
689 * and build start/count vectors. |
|
690 */ |
|
691 start[0] = tslice; |
|
692 count[0] = 1; |
|
693 for (i = 1; i < ndims; i++) { |
|
694 start[i] = 0; |
|
695 count[i] = (dim_idx[dimids[i]])->len; |
|
696 len *= (dim_idx[dimids[i]])->len; |
|
697 } |
|
698 |
|
699 val = alloc_var(nctype, len); |
|
700 |
|
701 switch (nctype) { |
|
702 case NC_FLOAT: |
|
703 wrap_nc(nc_get_vara_float(ncid, varid, start, count, val)); |
|
704 break; |
|
705 case NC_INT: |
|
706 wrap_nc(nc_get_vara_int(ncid, varid, start, count, val)); |
|
707 break; |
|
708 case NC_DOUBLE: |
|
709 wrap_nc(nc_get_vara_double(ncid, varid, start, count, val)); |
|
710 break; |
|
711 case NC_CHAR: |
|
712 wrap_nc(nc_get_vara_text(ncid, varid, start, count, val)); |
|
713 break; |
|
714 default: |
|
715 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
716 exit(3); |
|
717 } |
|
718 |
|
719 return val; |
|
720 } |
|
721 |
|
722 void write_var(int ncid, int varid, nc_type nctype, void *val) |
|
723 { |
|
724 switch (nctype) { |
|
725 case NC_FLOAT: |
|
726 wrap_nc(nc_put_var_float(ncid, varid, val)); |
|
727 break; |
|
728 case NC_INT: |
|
729 wrap_nc(nc_put_var_int(ncid, varid, val)); |
|
730 break; |
|
731 case NC_DOUBLE: |
|
732 wrap_nc(nc_put_var_double(ncid, varid, val)); |
|
733 break; |
|
734 case NC_CHAR: |
|
735 wrap_nc(nc_put_var_text(ncid, varid, val)); |
|
736 break; |
|
737 default: |
|
738 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
739 exit(3); |
|
740 } |
|
741 |
|
742 return; |
|
743 } |
|
744 |
|
745 void write_timeslice(int ncid, int varid, nc_type nctype, int ndims, int *dimids, struct dim **dim_idx, void *val, size_t tslice) |
|
746 { |
|
747 int i; |
|
748 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; |
|
749 |
|
750 /* Make sure time is really the first dimension */ |
|
751 if (strcmp((dim_idx[dimids[0]])->name, time_name)) { |
|
752 fprintf(stderr, "write_timeslice: %s is not first dimension of variable!\n", time_name); |
|
753 exit(3); |
|
754 } |
|
755 |
|
756 /* Build start/count vectors */ |
|
757 start[0] = tslice; |
|
758 count[0] = 1; |
|
759 for (i = 1; i < ndims; i++) { |
|
760 start[i] = 0; |
|
761 count[i] = (dim_idx[dimids[i]])->len; |
|
762 } |
|
763 |
|
764 switch (nctype) { |
|
765 case NC_FLOAT: |
|
766 wrap_nc(nc_put_vara_float(ncid, varid, start, count, val)); |
|
767 break; |
|
768 case NC_INT: |
|
769 wrap_nc(nc_put_vara_int(ncid, varid, start, count, val)); |
|
770 break; |
|
771 case NC_DOUBLE: |
|
772 wrap_nc(nc_put_vara_double(ncid, varid, start, count, val)); |
|
773 break; |
|
774 case NC_CHAR: |
|
775 wrap_nc(nc_put_vara_text(ncid, varid, start, count, val)); |
|
776 break; |
|
777 default: |
|
778 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
779 exit(3); |
|
780 } |
|
781 |
|
782 return; |
|
783 } |
|
784 |
|
785 void copy_metadata(int in_ncid, struct var *in_var_head, |
|
786 struct dim **in_dim_idx, int out_ncid, struct var *out_var_head) |
|
787 { |
|
788 int i, j; |
|
789 float hr[HOURS_PER_DAY]; |
|
790 struct var *in_vnode, *out_vnode; |
|
791 void *val; |
|
792 |
|
793 for (i = 0; i < HOURS_PER_DAY; i++) hr[i] = (float)i; |
|
794 |
|
795 if (in_var_head) { |
|
796 in_vnode = in_var_head; |
|
797 /* |
|
798 * March through input variables to stuff metadata values into |
|
799 * the new output files. NOTE: Time-based metadata variables |
|
800 * should contain only the last (time-slice) value (from all |
|
801 * input files). |
|
802 */ |
|
803 for (j = 0; j == 0 || in_vnode != in_var_head; j++) { |
|
804 if (in_vnode->metadata) { |
|
805 out_vnode = varlist_find_by_name(out_var_head, in_vnode->name); |
|
806 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { |
|
807 /* time is not the first dimension */ |
|
808 #ifdef DEBUG |
|
809 printf("Copying metadata variable %s\n", |
|
810 in_vnode->name); |
|
811 #endif /* DEBUG */ |
|
812 val = read_var(in_ncid, in_vnode->ncvarid, in_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx); |
|
813 write_var(out_ncid, out_vnode->ncvarid, out_vnode->nctype, val); |
|
814 } |
|
815 else { |
|
816 /* time is the first dimension */ |
|
817 #ifdef DEBUG |
|
818 printf("Copying last value of \ |
|
819 time-based metadata variable %s\n", in_vnode->name); |
|
820 #endif /* DEBUG */ |
|
821 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)); |
|
822 write_timeslice(out_ncid, out_vnode->ncvarid, out_vnode->nctype, in_vnode->ndims, in_vnode->dimids, in_dim_idx, val, 0); |
|
823 } |
|
824 free(val); |
|
825 /* |
|
826 * Just after the "time" variable, write out |
|
827 * the "hour" variable values. |
|
828 */ |
|
829 if (!strcmp(in_vnode->name, time_name)) { |
|
830 out_vnode = varlist_find_by_name(out_var_head, hour_name); |
|
831 write_var(out_ncid, out_vnode->ncvarid, |
|
832 out_vnode->nctype, hr); |
|
833 } |
|
834 } |
|
835 else { |
|
836 #ifdef DEBUG |
|
837 printf("Skipping variable %s\n", |
|
838 in_vnode->name); |
|
839 #endif /* DEBUG */ |
|
840 } |
|
841 in_vnode = in_vnode->next; |
|
842 } |
|
843 } |
|
844 |
|
845 return; |
|
846 } |
|
847 |
|
848 void open_inout(char *input_fname, char *mean_fname, char *stddev_fname, char *flist, size_t nsamples) |
|
849 { |
|
850 char *mean_history_gatt, *stddev_history_gatt, |
|
851 *mean_prefix="Statistical means from history files:", |
|
852 *stddev_prefix="Statistical standard deviations from history files:"; |
|
853 |
|
854 /* |
|
855 * Construct strings for history global attributes for the two output |
|
856 * files. |
|
857 */ |
|
858 if (!(mean_history_gatt = (char *)malloc((strlen(mean_prefix) + strlen(flist)+1)*sizeof(char)))) { |
|
859 perror("open_inout:mean_history_gatt"); |
|
860 exit(2); |
|
861 } |
|
862 if (!(stddev_history_gatt = (char *)malloc((strlen(stddev_prefix) + strlen(flist)+1)*sizeof(char)))) { |
|
863 perror("open_inout:stddev_history_gatt"); |
|
864 exit(2); |
|
865 } |
|
866 sprintf(mean_history_gatt, "%s%s", mean_prefix, flist); |
|
867 sprintf(stddev_history_gatt, "%s%s", stddev_prefix, flist); |
|
868 |
|
869 /* Open input file */ |
|
870 wrap_nc(nc_open(input_fname, NC_NOWRITE, &input_ncid)); |
|
871 /* Inquire about number of dimensions, variables, global attributes |
|
872 * and the ID of the unlimited dimension |
|
873 */ |
|
874 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, &input_ngatts, |
|
875 &input_unlimdimid)); |
|
876 /* Grab dimension IDs and lengths from input file */ |
|
877 get_dimensions(input_ncid, input_ndims, &input_dim_head, &input_dim_idx); |
|
878 /* Grab desired global attributes from input file */ |
|
879 get_gatts(input_ncid, input_ngatts, &input_gatt_head); |
|
880 /* Grab variable IDs and attributes from input file */ |
|
881 get_vars(input_ncid, input_nvars, &input_var_head); |
|
882 #ifdef DEBUG |
|
883 varlist_print(input_var_head); |
|
884 #endif /* DEBUG */ |
|
885 /* Create netCDF files for monthly means and standard deviations */ |
|
886 /* Create new netCDF files */ |
|
887 wrap_nc(nc_create(mean_fname, NC_NOCLOBBER, &mean_ncid)); |
|
888 wrap_nc(nc_create(stddev_fname, NC_NOCLOBBER, &stddev_ncid)); |
|
889 /* Define dimensions */ |
|
890 mean_ndims = put_dimensions(input_dim_head, input_ndims, |
|
891 input_unlimdimid, nsamples, &idid2mdid, mean_ncid, |
|
892 &mean_dim_head, &mean_unlimdimid, &mean_hour_dimid); |
|
893 stddev_ndims = put_dimensions(input_dim_head, input_ndims, |
|
894 input_unlimdimid, nsamples, &idid2sdid, stddev_ncid, |
|
895 &stddev_dim_head, &stddev_unlimdimid, &stddev_hour_dimid); |
|
896 /* Define variables and their attributes */ |
|
897 mean_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, |
|
898 idid2mdid, mean_ncid, mean_hour_dimid, &mean_var_head, |
|
899 MEAN_TYPE); |
|
900 stddev_nvars = define_vars(input_ncid, input_dim_idx, input_var_head, |
|
901 idid2sdid, stddev_ncid, stddev_hour_dimid, &stddev_var_head, |
|
902 STDDEV_TYPE); |
|
903 /* Store global attributes */ |
|
904 mean_ngatts = put_gatts(input_gatt_head, mean_ncid, mean_history_gatt); |
|
905 stddev_ngatts = put_gatts(input_gatt_head, stddev_ncid, |
|
906 stddev_history_gatt); |
|
907 /* End define mode */ |
|
908 wrap_nc(nc_enddef(mean_ncid)); |
|
909 wrap_nc(nc_enddef(stddev_ncid)); |
|
910 /* Write out metdata variables that do not depend on "time" */ |
|
911 copy_metadata(input_ncid, input_var_head, input_dim_idx, mean_ncid, |
|
912 mean_var_head); |
|
913 copy_metadata(input_ncid, input_var_head, input_dim_idx, stddev_ncid, |
|
914 stddev_var_head); |
|
915 |
|
916 wrap_nc(nc_close(input_ncid)); |
|
917 |
|
918 return; |
|
919 } |
|
920 |
|
921 size_t count_samples(int nifnames, char **input_fnames) |
|
922 { |
|
923 int i; |
|
924 char name[NC_MAX_NAME+1]; |
|
925 size_t len, cnt = 0; |
|
926 |
|
927 for (i = 0; i < nifnames; i++) { |
|
928 /* Open input file */ |
|
929 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); |
|
930 /* |
|
931 * Inquire about number of dimensions, variables, global |
|
932 * attributes and the ID of the unlimited dimension |
|
933 */ |
|
934 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, |
|
935 &input_ngatts, &input_unlimdimid)); |
|
936 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); |
|
937 if (strcmp(name, time_name)) { |
|
938 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); |
|
939 exit(4); |
|
940 } |
|
941 #ifdef DEBUG |
|
942 printf("%ld samples in %s\n", (long)len, input_fnames[i]); |
|
943 #endif /* DEBUG */ |
|
944 wrap_nc(nc_close(input_ncid)); |
|
945 cnt += len; |
|
946 } |
|
947 return cnt; |
|
948 } |
|
949 |
|
950 void alloc_means_stddevs(size_t d1, size_t d2, double ***meansp, double ***stddevsp, size_t ***cell_samplesp) |
|
951 { |
|
952 /* |
|
953 * Allocate space for arrays of means and standard deviations by |
|
954 * hour of day. |
|
955 */ |
|
956 int i; |
|
957 size_t **cell_samples; |
|
958 double **means, **stddevs; |
|
959 |
|
960 if (!(*meansp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { |
|
961 perror("alloc_means_stddevs:*meansp"); |
|
962 exit(2); |
|
963 } |
|
964 if (!(*stddevsp = (double **)malloc(HOURS_PER_DAY * sizeof(double *)))) { |
|
965 perror("alloc_means_stddevs:*stddevsp"); |
|
966 exit(2); |
|
967 } |
|
968 if (!(*cell_samplesp = (size_t **)malloc(HOURS_PER_DAY * sizeof(size_t *)))) { |
|
969 perror("alloc_means_stddevs:*cell_samplesp"); |
|
970 exit(2); |
|
971 } |
|
972 |
|
973 means = *meansp; |
|
974 stddevs = *stddevsp; |
|
975 cell_samples = *cell_samplesp; |
|
976 |
|
977 for (i = 0; i < HOURS_PER_DAY; i++) { |
|
978 if (!(means[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { |
|
979 perror("alloc_means_stddevs:means[i]"); |
|
980 exit(2); |
|
981 } |
|
982 if (!(stddevs[i] = (double *)malloc(d1 * d2 * sizeof(double)))) { |
|
983 perror("alloc_means_stddevs:stddevs[i]"); |
|
984 exit(2); |
|
985 } |
|
986 if (!(cell_samples[i] = (size_t *)malloc(d1 * d2 * sizeof(size_t)))) { |
|
987 perror("alloc_means_stddevs:cell_samples[i]"); |
|
988 exit(2); |
|
989 } |
|
990 } |
|
991 |
|
992 return; |
|
993 } |
|
994 |
|
995 void init_means_stddevs(size_t d1, size_t d2, double **means, double **stddevs, size_t **cell_samples, float FillValue) |
|
996 { |
|
997 int hr, i, j, pos; |
|
998 |
|
999 for (hr = 0; hr < HOURS_PER_DAY; hr++) { |
|
1000 for (i = 0; i < d1; i++) { |
|
1001 #pragma _CRI concurrent |
|
1002 for (j = 0; j < d2; j++) { |
|
1003 pos = i*d2+j; |
|
1004 means[hr][pos] = FillValue; |
|
1005 stddevs[hr][pos] = FillValue; |
|
1006 cell_samples[hr][pos] = 0; |
|
1007 } |
|
1008 } |
|
1009 } |
|
1010 |
|
1011 return; |
|
1012 } |
|
1013 |
|
1014 void reset_cell_samples(size_t d1, size_t d2, size_t **cell_samples) |
|
1015 { |
|
1016 int i, j, hr, pos; |
|
1017 |
|
1018 for (hr = 0; hr < HOURS_PER_DAY; hr++) { |
|
1019 for (i = 0; i < d1; i++) { |
|
1020 #pragma _CRI concurrent |
|
1021 for (j = 0; j < d2; j++) { |
|
1022 pos = i*d2+j; |
|
1023 cell_samples[hr][pos] = 0; |
|
1024 } |
|
1025 } |
|
1026 } |
|
1027 |
|
1028 return; |
|
1029 } |
|
1030 |
|
1031 void free_means_stddevs(double **means, double **stddevs, size_t **cell_samples) |
|
1032 { |
|
1033 /* |
|
1034 * Free space from arrays of means and standard deviations by |
|
1035 * hour of day. |
|
1036 */ |
|
1037 int i; |
|
1038 |
|
1039 for (i = 0; i < HOURS_PER_DAY; i++) { |
|
1040 free(means[i]); |
|
1041 free(stddevs[i]); |
|
1042 free(cell_samples[i]); |
|
1043 } |
|
1044 |
|
1045 free(means); |
|
1046 free(stddevs); |
|
1047 free(cell_samples); |
|
1048 |
|
1049 return; |
|
1050 } |
|
1051 |
|
1052 void add_to_means_float(float *val, int sec, size_t d1, size_t d2, |
|
1053 char FillFlag, float FillValue, double **means, size_t **cell_samples) |
|
1054 { |
|
1055 int i, j, hr, pos; |
|
1056 |
|
1057 hr = (int)floor((double)sec/(double)SEC_PER_HOUR); |
|
1058 |
|
1059 for (i = 0; i < d1; i++) { |
|
1060 #pragma _CRI concurrent |
|
1061 for (j = 0; j < d2; j++) { |
|
1062 pos = i*d2+j; |
|
1063 if (!FillFlag || (FillFlag && val[pos] != FillValue)) { |
|
1064 if (cell_samples[hr][pos] == 0) |
|
1065 means[hr][pos] = (double)val[pos]; |
|
1066 else |
|
1067 means[hr][pos] += (double)val[pos]; |
|
1068 ++cell_samples[hr][pos]; |
|
1069 } |
|
1070 } |
|
1071 } |
|
1072 |
|
1073 return; |
|
1074 } |
|
1075 |
|
1076 void add_to_means_double(double *val, int sec, size_t d1, size_t d2, |
|
1077 char FillFlag, float FillValue, double **means, size_t **cell_samples) |
|
1078 { |
|
1079 int i, j, hr, pos; |
|
1080 |
|
1081 hr = (int)floor((double)sec/(double)SEC_PER_HOUR); |
|
1082 |
|
1083 for (i = 0; i < d1; i++) { |
|
1084 #pragma _CRI concurrent |
|
1085 for (j = 0; j < d2; j++) { |
|
1086 pos = i*d2+j; |
|
1087 if (!FillFlag || (FillFlag && val[pos] != FillValue)) { |
|
1088 if (cell_samples[hr][pos] == 0) |
|
1089 means[hr][pos] = val[pos]; |
|
1090 else |
|
1091 means[hr][pos] += val[pos]; |
|
1092 ++cell_samples[hr][pos]; |
|
1093 } |
|
1094 } |
|
1095 } |
|
1096 |
|
1097 return; |
|
1098 } |
|
1099 |
|
1100 |
|
1101 void divide_means(size_t d1, size_t d2, double **means, size_t **cell_samples) |
|
1102 { |
|
1103 int i, j, hr, pos; |
|
1104 |
|
1105 for (hr = 0; hr < HOURS_PER_DAY; hr++) { |
|
1106 for (i = 0; i < d1; i++) { |
|
1107 #pragma _CRI concurrent |
|
1108 for (j = 0; j < d2; j++) { |
|
1109 pos = i*d2+j; |
|
1110 if (cell_samples[hr][pos] != 0) { |
|
1111 means[hr][pos] /= (double)cell_samples[hr][pos]; |
|
1112 } |
|
1113 } |
|
1114 } |
|
1115 } |
|
1116 |
|
1117 return; |
|
1118 } |
|
1119 |
|
1120 void add_to_stddevs_float(float *val, int sec, size_t d1, size_t d2, |
|
1121 char FillFlag, float FillValue, double **means, double **stddevs, |
|
1122 size_t **cell_samples) |
|
1123 { |
|
1124 int i, j, hr, pos; |
|
1125 double delta; |
|
1126 |
|
1127 hr = (int)floor((double)sec/(double)SEC_PER_HOUR); |
|
1128 |
|
1129 for (i = 0; i < d1; i++) { |
|
1130 #pragma _CRI concurrent |
|
1131 for (j = 0; j < d2; j++) { |
|
1132 pos = i*d2+j; |
|
1133 if (!FillFlag || (FillFlag && val[pos] != FillValue |
|
1134 && means[hr][pos] != FillValue)) { |
|
1135 delta = means[hr][pos] - (double)val[pos]; |
|
1136 if (cell_samples[hr][pos] == 0) |
|
1137 stddevs[hr][pos] = delta * delta; |
|
1138 else |
|
1139 stddevs[hr][pos] += delta * delta; |
|
1140 ++cell_samples[hr][pos]; |
|
1141 } |
|
1142 } |
|
1143 } |
|
1144 |
|
1145 return; |
|
1146 } |
|
1147 |
|
1148 void add_to_stddevs_double(double *val, int sec, size_t d1, size_t d2, |
|
1149 char FillFlag, float FillValue, double **means, double **stddevs, |
|
1150 size_t **cell_samples) |
|
1151 { |
|
1152 int i, j, hr, pos; |
|
1153 double delta; |
|
1154 |
|
1155 hr = (int)floor((double)sec/(double)SEC_PER_HOUR); |
|
1156 |
|
1157 for (i = 0; i < d1; i++) { |
|
1158 #pragma _CRI concurrent |
|
1159 for (j = 0; j < d2; j++) { |
|
1160 pos = i*d2+j; |
|
1161 if (!FillFlag || (FillFlag && val[pos] != FillValue |
|
1162 && means[hr][pos] != FillValue)) { |
|
1163 delta = means[hr][pos] - val[pos]; |
|
1164 if (cell_samples[hr][pos] == 0) |
|
1165 stddevs[hr][pos] = delta * delta; |
|
1166 else |
|
1167 stddevs[hr][pos] += delta * delta; |
|
1168 ++cell_samples[hr][pos]; |
|
1169 } |
|
1170 } |
|
1171 } |
|
1172 |
|
1173 return; |
|
1174 } |
|
1175 |
|
1176 |
|
1177 void divide_sqrt_stddevs(size_t d1, size_t d2, double **stddevs, size_t **cell_samples) |
|
1178 { |
|
1179 int i, j, hr, pos; |
|
1180 |
|
1181 for (hr = 0; hr < HOURS_PER_DAY; hr++) { |
|
1182 for (i = 0; i < d1; i++) { |
|
1183 #pragma _CRI concurrent |
|
1184 for (j = 0; j < d2; j++) { |
|
1185 pos = i*d2+j; |
|
1186 if (cell_samples[hr][pos] != 0) { |
|
1187 stddevs[hr][pos] /= (double)cell_samples[hr][pos]; |
|
1188 stddevs[hr][pos] = sqrt(stddevs[hr][pos]); |
|
1189 } |
|
1190 } |
|
1191 } |
|
1192 } |
|
1193 |
|
1194 return; |
|
1195 } |
|
1196 |
|
1197 void alloc_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void ***valsp, int **mcsecp) |
|
1198 { |
|
1199 int i; |
|
1200 void **vals; |
|
1201 |
|
1202 /* Allocate space for all timeslices */ |
|
1203 if (!(*valsp = (void **)malloc(nsamples * sizeof(void *)))) { |
|
1204 perror("alloc_all_samples:*valsp"); |
|
1205 exit(2); |
|
1206 } |
|
1207 vals = *valsp; |
|
1208 for (i = 0; i < nsamples; i++) { |
|
1209 vals[i] = alloc_var(nctype, (d1 * d2)); |
|
1210 } |
|
1211 if (!(*mcsecp = (int *)malloc(nsamples * sizeof(int)))) { |
|
1212 perror("alloc_all_samples:*mcsecp"); |
|
1213 exit(2); |
|
1214 } |
|
1215 |
|
1216 return; |
|
1217 } |
|
1218 |
|
1219 void init_all_samples(size_t d1, size_t d2, size_t nsamples, nc_type nctype, void **vals, int *mcsec) |
|
1220 { |
|
1221 int i, j; |
|
1222 float **fvals = NULL; |
|
1223 double **dvals = NULL; |
|
1224 |
|
1225 switch(nctype) { |
|
1226 case NC_FLOAT: |
|
1227 fvals = (float **)vals; |
|
1228 break; |
|
1229 case NC_DOUBLE: |
|
1230 dvals = (double **)vals; |
|
1231 break; |
|
1232 default: |
|
1233 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
1234 exit(3); |
|
1235 } |
|
1236 |
|
1237 for (i = 0; i < nsamples; i++) { |
|
1238 for (j = 0; j < (d1 * d2); j++) { |
|
1239 switch(nctype) { |
|
1240 case NC_FLOAT: |
|
1241 fvals[i][j] = (float)0.0; |
|
1242 break; |
|
1243 case NC_DOUBLE: |
|
1244 dvals[i][j] = (double)0.0; |
|
1245 break; |
|
1246 default: |
|
1247 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
1248 exit(3); |
|
1249 } |
|
1250 } |
|
1251 mcsec[i] = 0; |
|
1252 } |
|
1253 |
|
1254 return; |
|
1255 } |
|
1256 |
|
1257 void free_all_samples(size_t nsamples, void **vals, int *mcsec) |
|
1258 { |
|
1259 int i; |
|
1260 |
|
1261 for (i = 0; i < nsamples; i++) |
|
1262 if (vals[i]) free(vals[i]); |
|
1263 |
|
1264 free(vals); |
|
1265 free(mcsec); |
|
1266 |
|
1267 return; |
|
1268 } |
|
1269 |
|
1270 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) |
|
1271 { |
|
1272 int i, ts; |
|
1273 int varid; |
|
1274 size_t len, nslice = 0, start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; |
|
1275 char name[NC_MAX_NAME+1]; |
|
1276 |
|
1277 for (i = 0; i < nifnames; i++) { |
|
1278 #ifdef DEBUG |
|
1279 printf("\tOpening %s", input_fnames[i]); |
|
1280 if (ndims > 3) printf(" to retrieve level %d\n", level); |
|
1281 #endif /* DEBUG */ |
|
1282 /* Open input file */ |
|
1283 wrap_nc(nc_open(input_fnames[i], NC_NOWRITE, &input_ncid)); |
|
1284 /* |
|
1285 * Inquire about number of dimensions, variables, global |
|
1286 * attributes and the ID of the unlimited dimension |
|
1287 */ |
|
1288 wrap_nc(nc_inq(input_ncid, &input_ndims, &input_nvars, |
|
1289 &input_ngatts, &input_unlimdimid)); |
|
1290 wrap_nc(nc_inq_dim(input_ncid, input_unlimdimid, name, &len)); |
|
1291 if (strcmp(name, time_name)) { |
|
1292 fprintf(stderr, "%s is not the unlimited dimension for file %s!\n", time_name, input_fnames[i]); |
|
1293 exit(4); |
|
1294 } |
|
1295 /* Make sure we don't run off the end of allocated space */ |
|
1296 if ((nslice+len) > nsamples) { |
|
1297 fprintf(stderr, "Number of time slices exceeds previous count of %ld!\n", (long)nsamples); |
|
1298 exit(4); |
|
1299 } |
|
1300 /* Get seconds of day variable */ |
|
1301 wrap_nc(nc_inq_varid(input_ncid, mcsec_name, &varid)); |
|
1302 wrap_nc(nc_get_var_int(input_ncid, varid, &mcsec[nslice])); |
|
1303 /* Get variable ID for requested variable */ |
|
1304 wrap_nc(nc_inq_varid(input_ncid, var_name, &varid)); |
|
1305 /* Retrieve time slice of desired variable */ |
|
1306 for (ts = 0; ts < len; ts++) { |
|
1307 start[0] = ts; |
|
1308 count[0] = 1; |
|
1309 start[(ndims-2)] = 0; |
|
1310 count[(ndims-2)] = d1; |
|
1311 start[(ndims-1)] = 0; |
|
1312 count[(ndims-1)] = d2; |
|
1313 if (ndims == 4) { |
|
1314 start[1] = level; |
|
1315 count[1] = 1; |
|
1316 } |
|
1317 switch(nctype) { |
|
1318 case NC_FLOAT: |
|
1319 wrap_nc(nc_get_vara_float(input_ncid, varid, start, count, vals[nslice])); |
|
1320 break; |
|
1321 case NC_DOUBLE: |
|
1322 wrap_nc(nc_get_vara_double(input_ncid, varid, start, count, vals[nslice])); |
|
1323 break; |
|
1324 default: |
|
1325 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
1326 exit(3); |
|
1327 } |
|
1328 nslice++; |
|
1329 } |
|
1330 |
|
1331 /* Close input file */ |
|
1332 wrap_nc(nc_close(input_ncid)); |
|
1333 } |
|
1334 |
|
1335 return; |
|
1336 } |
|
1337 |
|
1338 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) |
|
1339 { |
|
1340 int i; |
|
1341 |
|
1342 for (i = 0; i < nsamples; i++) { |
|
1343 switch(nctype) { |
|
1344 case NC_FLOAT: |
|
1345 if (stat_type == MEAN_TYPE) |
|
1346 add_to_means_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples); |
|
1347 else |
|
1348 add_to_stddevs_float(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); |
|
1349 break; |
|
1350 case NC_DOUBLE: |
|
1351 if (stat_type == MEAN_TYPE) |
|
1352 add_to_means_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, cell_samples); |
|
1353 else |
|
1354 add_to_stddevs_double(vals[i], mcsec[i], d1, d2, FillFlag, FillValue, means, stddevs, cell_samples); |
|
1355 break; |
|
1356 default: |
|
1357 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
1358 exit(3); |
|
1359 } |
|
1360 } |
|
1361 |
|
1362 /* Divide sums by number of samples */ |
|
1363 if (stat_type == MEAN_TYPE) |
|
1364 divide_means(d1, d2, means, cell_samples); |
|
1365 else |
|
1366 divide_sqrt_stddevs(d1, d2, stddevs, cell_samples); |
|
1367 |
|
1368 return; |
|
1369 } |
|
1370 |
|
1371 float *double_to_float(size_t len, double *dvar) |
|
1372 { |
|
1373 int i; |
|
1374 float *fvar; |
|
1375 |
|
1376 if (!(fvar = (float *)malloc(len * sizeof(float)))) { |
|
1377 perror("double_to_float:fvar"); |
|
1378 exit(2); |
|
1379 } |
|
1380 |
|
1381 for (i = 0; i < len; i++) |
|
1382 fvar[i] = (float)dvar[i]; |
|
1383 |
|
1384 return fvar; |
|
1385 } |
|
1386 |
|
1387 void write_var_hours(int ncid, int varid, nc_type nctype, size_t d1, size_t d2, |
|
1388 int level, int ndims, double **var_hours) |
|
1389 { |
|
1390 /* Output dimensions should be one larger than input dimensions */ |
|
1391 int i, hr; |
|
1392 size_t start[NC_MAX_VAR_DIMS], count[NC_MAX_VAR_DIMS]; |
|
1393 float *fvar = NULL; |
|
1394 |
|
1395 if (nctype == NC_FLOAT) { |
|
1396 if (!(fvar = (float *)malloc(d1 * d2 * sizeof(float)))) { |
|
1397 perror("write_var_hours:fvar"); |
|
1398 exit(2); |
|
1399 } |
|
1400 } |
|
1401 |
|
1402 for (hr = 0; hr < HOURS_PER_DAY; hr++) { |
|
1403 start[0] = 0; |
|
1404 count[0] = 1; /* time */ |
|
1405 start[1] = hr; |
|
1406 count[1] = 1; /* hour */ |
|
1407 start[(ndims-2)] = 0; |
|
1408 count[(ndims-2)] = d1; |
|
1409 start[(ndims-1)] = 0; |
|
1410 count[(ndims-1)] = d2; |
|
1411 if (ndims == 5) { |
|
1412 start[2] = level; |
|
1413 count[2] = 1; |
|
1414 } |
|
1415 switch (nctype) { |
|
1416 case NC_FLOAT: |
|
1417 for (i = 0; i < (d1 * d2); i++) |
|
1418 fvar[i] = (float)var_hours[hr][i]; |
|
1419 wrap_nc(nc_put_vara_float(ncid, varid, start, count, fvar)); |
|
1420 break; |
|
1421 case NC_DOUBLE: |
|
1422 wrap_nc(nc_put_vara_double(ncid, varid, start, count, var_hours[hr])); |
|
1423 break; |
|
1424 default: |
|
1425 fprintf(stderr, "netCDF external data type %d not supported\n", nctype); |
|
1426 exit(3); |
|
1427 } |
|
1428 } |
|
1429 |
|
1430 if (nctype == NC_FLOAT) |
|
1431 free(fvar); |
|
1432 |
|
1433 return; |
|
1434 } |
|
1435 |
|
1436 void compute_stats(int nifnames, char **input_fnames, size_t nsamples) |
|
1437 { |
|
1438 int j, k, nlevels, *mcsec; |
|
1439 size_t d1, d2, **cell_samples; |
|
1440 double **means, **stddevs; |
|
1441 struct var *in_vnode, *out_vnode; |
|
1442 void **vals; |
|
1443 |
|
1444 if (input_var_head) { |
|
1445 in_vnode = input_var_head; |
|
1446 /* March through non-metadata variables to compute stats */ |
|
1447 for (j = 0; j == 0 || in_vnode != input_var_head; j++) { |
|
1448 if (!in_vnode->metadata) { |
|
1449 /* Make sure time is really the first dimension */ |
|
1450 if (strcmp((input_dim_idx[in_vnode->dimids[0]])->name, time_name)) { |
|
1451 fprintf(stderr, "compute_stats: %s is not first dimension of variable %s!\n", time_name, in_vnode->name); |
|
1452 exit(3); |
|
1453 } |
|
1454 /* Figure out input dimensions */ |
|
1455 if (in_vnode->ndims < 3 || in_vnode->ndims > 4) { |
|
1456 fprintf(stderr, "compute_stats: %s has %d dimensions!\n", in_vnode->name, in_vnode->ndims); |
|
1457 exit(3); |
|
1458 } |
|
1459 /* Assume 2D output; find dimensions */ |
|
1460 d1 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-2)]])->len; |
|
1461 d2 = (input_dim_idx[in_vnode->dimids[(in_vnode->ndims-1)]])->len; |
|
1462 if (in_vnode->ndims == 3) |
|
1463 nlevels = 1; |
|
1464 else |
|
1465 nlevels = (input_dim_idx[in_vnode->dimids[1]])->len; |
|
1466 /* Allocate working space for means and stddevs */ |
|
1467 alloc_means_stddevs(d1, d2, &means, &stddevs, &cell_samples); |
|
1468 init_means_stddevs(d1, d2, means, stddevs, cell_samples, in_vnode->FillValue); |
|
1469 /* Allocate and initize space for entire field across all files */ |
|
1470 alloc_all_samples(d1, d2, nsamples, in_vnode->nctype, &vals, &mcsec); |
|
1471 init_all_samples(d1, d2, nsamples, in_vnode->nctype, vals, mcsec); |
|
1472 printf("Computing statistics for %s\n", |
|
1473 in_vnode->name); |
|
1474 /* Compute means and stddevs, then write them */ |
|
1475 for (k = 0; k < nlevels; k++) { |
|
1476 /* Read the entire fields from all files */ |
|
1477 read_all_samples(nifnames, input_fnames, d1, d2, nsamples, in_vnode->name, in_vnode->nctype, k, in_vnode->ndims, vals, mcsec); |
|
1478 /* Compute means */ |
|
1479 compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, MEAN_TYPE, means, stddevs, cell_samples); |
|
1480 /* Find corresponding output variable on the mean output file */ |
|
1481 out_vnode = varlist_find_by_name(mean_var_head, in_vnode->name); |
|
1482 /* Write out the means for this variable */ |
|
1483 write_var_hours(mean_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, means); |
|
1484 /* Zero out cell_samples so they can be used as a flag again for computing stddevs */ |
|
1485 reset_cell_samples(d1, d2, cell_samples); |
|
1486 /* Compute stddevs using means */ |
|
1487 compute(d1, d2, nsamples, in_vnode->nctype, vals, mcsec, in_vnode->FillFlag, in_vnode->FillValue, STDDEV_TYPE, means, stddevs, cell_samples); |
|
1488 /* Find corresponding output variable on the stddev output file */ |
|
1489 out_vnode = varlist_find_by_name(stddev_var_head, in_vnode->name); |
|
1490 /* Write out stddevs for this variable */ |
|
1491 write_var_hours(stddev_ncid, out_vnode->ncvarid, out_vnode->nctype, d1, d2, k, out_vnode->ndims, stddevs); |
|
1492 } |
|
1493 |
|
1494 /* Free all samples */ |
|
1495 free_all_samples(nsamples, vals, mcsec); |
|
1496 /* Free working space for means and stddevs */ |
|
1497 free_means_stddevs(means, stddevs, cell_samples); |
|
1498 } |
|
1499 in_vnode = in_vnode->next; |
|
1500 } |
|
1501 } |
|
1502 return; |
|
1503 } |
|
1504 |
|
1505 void usage(char *program) |
|
1506 { |
|
1507 fprintf(stderr, "Usage: %s -m mean.nc -s stddev.nc hist_file1.nc [ hist_file2.nc ... ]\n", program); |
|
1508 fprintf(stderr, "WARNING: It is assumed that the list of input files is specified in time order!\n"); |
|
1509 exit(4); |
|
1510 } |
|
1511 |
|
1512 int main(int argc, char **argv) |
|
1513 { |
|
1514 int i, ic, nifnames; |
|
1515 size_t len, pos, nsamples; |
|
1516 char *mfname, *sfname, **ifnames, *flist; |
|
1517 |
|
1518 ifnames = NULL; |
|
1519 mfname = sfname = flist = NULL; |
|
1520 nifnames = 0; |
|
1521 |
|
1522 #ifdef DEBUG |
|
1523 printf("Number of metadata variables (nmvars) = %d\n", nmvars); |
|
1524 fflush(stdout); |
|
1525 #endif /* DEBUG */ |
|
1526 |
|
1527 /* check command-line flags and switches */ |
|
1528 while ((ic = getopt(argc, argv, "m:s:")) != -1) { |
|
1529 switch(ic) { |
|
1530 case 'm': |
|
1531 if (!(mfname = strdup(optarg))) { |
|
1532 perror("mfname"); |
|
1533 exit(2); |
|
1534 } |
|
1535 break; |
|
1536 case 's': |
|
1537 if (!(sfname = strdup(optarg))) { |
|
1538 perror("sfname"); |
|
1539 exit(2); |
|
1540 } |
|
1541 break; |
|
1542 } |
|
1543 } |
|
1544 |
|
1545 if (!mfname) { |
|
1546 fprintf(stderr, "Output file name for writing means required.\n"); |
|
1547 usage(argv[0]); |
|
1548 } |
|
1549 if (!sfname) { |
|
1550 fprintf(stderr, "Output file name for writing standard deviations required.\n"); |
|
1551 usage(argv[0]); |
|
1552 } |
|
1553 |
|
1554 if (optind < argc) { |
|
1555 if (!(ifnames = (char **)malloc((argc-optind+1)*sizeof(char *)))) { |
|
1556 perror("ifnames"); |
|
1557 exit(2); |
|
1558 } |
|
1559 for (i = optind; i < argc; i++) { |
|
1560 if (!(ifnames[nifnames++] = strdup(argv[i]))) { |
|
1561 perror("ifnames[nifnames++]"); |
|
1562 exit(2); |
|
1563 } |
|
1564 } |
|
1565 } |
|
1566 else { |
|
1567 fprintf(stderr, "At least one input file name is required.\n"); |
|
1568 usage(argv[0]); |
|
1569 } |
|
1570 |
|
1571 |
|
1572 /* |
|
1573 * Build list of input files to be included in the global history |
|
1574 * attribute on the two outputs files. |
|
1575 */ |
|
1576 for (i = len = 0; i < nifnames; i++) |
|
1577 len += strlen(ifnames[i]); |
|
1578 len += nifnames + 1; /* space in front of every name + null terminator */ |
|
1579 if (!(flist = (char *)malloc(len * sizeof(char)))) { |
|
1580 perror("flist"); |
|
1581 exit(2); |
|
1582 } |
|
1583 for (i = pos = 0; i < nifnames; pos += (strlen(ifnames[i])+1), i++) |
|
1584 sprintf(&flist[pos], " %s", ifnames[i]); |
|
1585 #ifdef DEBUG |
|
1586 printf("flist=%s\n", flist); |
|
1587 fflush(stdout); |
|
1588 #endif /* DEBUG */ |
|
1589 |
|
1590 /* Open every file to count up the number of time samples. */ |
|
1591 nsamples = count_samples(nifnames, ifnames); |
|
1592 #ifdef DEBUG |
|
1593 printf("Number of samples across all files = %ld\n", (long)nsamples); |
|
1594 fflush(stdout); |
|
1595 #endif /* DEBUG */ |
|
1596 |
|
1597 /* |
|
1598 * Open the *last* input file and the two output files (for mean and |
|
1599 * standard deviation). Define dimensions, variables, and attributes |
|
1600 * for the two output files based on the *last* input files. The |
|
1601 * last file is used because some metadata variables will contain |
|
1602 * only values from the last time slice from the period over which |
|
1603 * the statistics are computed. |
|
1604 */ |
|
1605 open_inout(ifnames[(nifnames-1)], mfname, sfname, flist, nsamples); |
|
1606 |
|
1607 compute_stats(nifnames, ifnames, nsamples); |
|
1608 |
|
1609 /* Close the two output files */ |
|
1610 wrap_nc(nc_close(mean_ncid)); |
|
1611 wrap_nc(nc_close(stddev_ncid)); |
|
1612 |
|
1613 return 0; |
|
1614 } |