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