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