1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 CMORlite cmorlite CMOR lite
12 */
13
14 #include <cdi.h>
15
16 #include "cdo_options.h"
17 #include "process_int.h"
18 #include "param_conversion.h"
19 #include "cdo_cmor.h"
20 #include "pmlist.h"
21 #include "convert_units.h"
22 #include "cdo_cdi_wrapper.h"
23 #include "cdi_lockedIO.h"
24
25 void
cdo_define_var_units(CmorVar & var,const int vlistID2,const int varID,const char * const units)26 cdo_define_var_units(CmorVar &var, const int vlistID2, const int varID, const char *const units)
27 {
28 char units_old[CDI_MAX_NAME];
29
30 vlistInqVarUnits(vlistID2, varID, units_old);
31 const auto len1 = strlen(units_old);
32 const auto len2 = strlen(units);
33
34 if (strcmp(units, units_old) != 0)
35 {
36 if (len1 > 0 && len2 > 0)
37 {
38 var.changeunits = true;
39 strcpy(var.units_old, units_old);
40 strcpy(var.units, units);
41 }
42
43 cdiDefKeyString(vlistID2, varID, CDI_KEY_UNITS, units);
44 cdiDefAttTxt(vlistID2, varID, "original_units", (int) strlen(units_old), units_old);
45 }
46 }
47
48 void
cmor_check_init(int nvars,std::vector<CmorVar> & vars)49 cmor_check_init(int nvars, std::vector<CmorVar> &vars)
50 {
51 for (int varID = 0; varID < nvars; ++varID)
52 {
53 auto &var = vars[varID];
54 if (var.checkvalid || var.check_min_mean_abs || var.check_max_mean_abs)
55 {
56 var.amean = 0;
57 var.nvals = 0;
58 var.n_lower_min = 0;
59 var.n_greater_max = 0;
60 }
61 }
62 }
63
64 void
cmor_check_eval(int vlistID,int nvars,const std::vector<CmorVar> & vars)65 cmor_check_eval(int vlistID, int nvars, const std::vector<CmorVar> &vars)
66 {
67 char varname[CDI_MAX_NAME];
68
69 for (int varID = 0; varID < nvars; ++varID)
70 {
71 const auto &var = vars[varID];
72 if (var.checkvalid || var.check_min_mean_abs || var.check_max_mean_abs)
73 {
74 auto amean = var.amean;
75 const auto nvals = var.nvals;
76
77 if (nvals > 0) amean /= nvals;
78
79 const auto n_lower_min = var.n_lower_min;
80 const auto n_greater_max = var.n_greater_max;
81
82 vlistInqVarName(vlistID, varID, varname);
83
84 if (n_lower_min > 0)
85 cdo_warning("Invalid value(s) detected for variable '%s': %ld values were lower than minimum valid value (%.4g).",
86 varname, n_lower_min, var.valid_min);
87 if (n_greater_max > 0)
88 cdo_warning("Invalid value(s) detected for variable '%s': %ld values were greater than maximum valid value (%.4g).",
89 varname, n_greater_max, var.valid_max);
90
91 if (var.check_min_mean_abs)
92 {
93 if (amean < .1 * var.ok_min_mean_abs)
94 cdo_warning("Invalid Absolute Mean for variable '%s' (%.5g) is lower by more than an order of magnitude than "
95 "minimum allowed: %.4g",
96 varname, amean, var.ok_min_mean_abs);
97
98 if (amean < var.ok_min_mean_abs)
99 cdo_warning("Invalid Absolute Mean for variable '%s' (%.5g) is lower than minimum allowed: %.4g", varname, amean,
100 var.ok_min_mean_abs);
101 }
102
103 if (var.check_max_mean_abs)
104 {
105 if (amean > 10. * var.ok_max_mean_abs)
106 cdo_warning("Invalid Absolute Mean for variable '%s' (%.5g) is greater by more than an order of magnitude than "
107 "maximum allowed: %.4g",
108 varname, amean, var.ok_max_mean_abs);
109
110 if (amean > var.ok_max_mean_abs)
111 cdo_warning("Invalid Absolute Mean for variable '%s' (%.5g) is greater than maximum allowed: %.4g", varname, amean,
112 var.ok_max_mean_abs);
113 }
114 }
115 }
116 }
117
118 void
cmor_check_prep(CmorVar & var,const long gridsize,const double missval,const double * const array)119 cmor_check_prep(CmorVar &var, const long gridsize, const double missval, const double *const array)
120 {
121 if (var.checkvalid || var.check_min_mean_abs || var.check_max_mean_abs)
122 {
123 double amean = 0;
124 long nvals = 0;
125
126 for (long i = 0; i < gridsize; ++i)
127 {
128 const auto aval = array[i];
129 if (!DBL_IS_EQUAL(aval, missval))
130 {
131 amean += std::fabs(aval);
132 nvals++;
133 }
134 }
135
136 var.amean += amean;
137 var.nvals += nvals;
138
139 long n_lower_min = 0, n_greater_max = 0;
140
141 for (long i = 0; i < gridsize; ++i)
142 {
143 const auto aval = array[i];
144 if (!DBL_IS_EQUAL(aval, missval))
145 {
146 if (aval < var.valid_min) n_lower_min++;
147 if (aval > var.valid_max) n_greater_max++;
148 }
149 }
150
151 var.n_lower_min += n_lower_min;
152 var.n_greater_max += n_greater_max;
153 }
154 }
155
156 static void
applyCmorList(PMList & pmlist,int nvars,int vlistID2,std::vector<CmorVar> & vars)157 applyCmorList(PMList &pmlist, int nvars, int vlistID2, std::vector<CmorVar> &vars)
158 {
159 const std::vector<std::string> hentry = { "Header" };
160 const std::vector<std::string> ventry = { "variable_entry", "parameter" };
161 char varname[CDI_MAX_NAME];
162
163 // search for global missing value
164 bool lmissval = false;
165 double missval = 0;
166
167 {
168 auto kvlist = pmlist.getKVListVentry(hentry);
169 if (kvlist)
170 {
171 for (const auto &kv : *kvlist)
172 {
173 const auto &key = kv.key;
174 const auto &value = kv.values[0];
175 if (kv.nvalues != 1 || value.empty()) continue;
176
177 if (key == "missing_value")
178 {
179 lmissval = true;
180 missval = parameter_to_double(value);
181 }
182 else if (key == "table_id" || key == "modeling_realm" || key == "realm" || key == "project_id" || key == "frequency")
183 {
184 cdiDefAttTxt(vlistID2, CDI_GLOBAL, key.c_str(), (int) value.size(), value.c_str());
185 }
186 }
187 }
188 }
189
190 for (int varID = 0; varID < nvars; varID++)
191 {
192 auto &var = vars[varID];
193 vlistInqVarName(vlistID2, varID, varname);
194
195 strcpy(var.name, varname);
196 if (lmissval)
197 {
198 const auto missval_old = vlistInqVarMissval(vlistID2, varID);
199 if (!DBL_IS_EQUAL(missval, missval_old))
200 {
201 var.changemissval = true;
202 var.missval_old = missval_old;
203 vlistDefVarMissval(vlistID2, varID, missval);
204 }
205 }
206
207 auto kvlist = pmlist.searchKVListVentry("name", varname, ventry);
208 if (kvlist)
209 {
210 bool lvalid_min = false, lvalid_max = false;
211
212 for (const auto &kv : *kvlist)
213 {
214 const auto &key = kv.key;
215 const auto &value = kv.values[0];
216 if (kv.nvalues != 1 || value.empty()) continue;
217 auto value_cstr = value.c_str();
218
219 // printf("key=%s value=>%s<\n", key.c_str(), value.c_str());
220
221 // clang-format off
222 if (key == "standard_name") cdiDefKeyString(vlistID2, varID, CDI_KEY_STDNAME, value_cstr);
223 else if (key == "long_name") cdiDefKeyString(vlistID2, varID, CDI_KEY_LONGNAME, value_cstr);
224 else if (key == "units") cdo_define_var_units(var, vlistID2, varID, value_cstr);
225 else if (key == "name") ; // cdiDefKeyString(vlistID2, varID, CDI_KEY_NAME, parameter_to_word(value.c_str()));
226 else if (key == "out_name")
227 {
228 const auto outname = parameter_to_word(value.c_str());
229 if (!cdo_cmpstr(var.name, outname))
230 {
231 cdiDefKeyString(vlistID2, varID, CDI_KEY_NAME, outname);
232 cdiDefAttTxt(vlistID2, varID, "original_name", (int) strlen(var.name), var.name);
233 }
234 }
235 else if (key == "param") vlistDefVarParam(vlistID2, varID, string_to_param(parameter_to_word(value.c_str())));
236 else if (key == "out_param") vlistDefVarParam(vlistID2, varID, string_to_param(parameter_to_word(value.c_str())));
237 else if (key == "comment") cdiDefAttTxt(vlistID2, varID, key.c_str(), (int) value.size(), value_cstr);
238 else if (key == "cell_methods") cdiDefAttTxt(vlistID2, varID, key.c_str(), (int) value.size(), value_cstr);
239 else if (key == "cell_measures") cdiDefAttTxt(vlistID2, varID, key.c_str(), (int) value.size(), value_cstr);
240 else if (key == "delete") var.remove = parameter_to_bool(value);
241 else if (key == "convert") var.convert = parameter_to_bool(value);
242 else if (key == "factor")
243 {
244 var.lfactor = true;
245 var.factor = parameter_to_double(value);
246 if (Options::cdoVerbose) cdo_print("%s - scale factor %g", varname, var.factor);
247 }
248 else if (key == "missval" || key == "missing_value")
249 {
250 missval = parameter_to_double(value);
251 auto missval_old = vlistInqVarMissval(vlistID2, varID);
252 if (!DBL_IS_EQUAL(missval, missval_old))
253 {
254 if (Options::cdoVerbose) cdo_print("%s - change missval from %g to %g", varname, missval_old, missval);
255 var.changemissval = true;
256 var.missval_old = missval_old;
257 vlistDefVarMissval(vlistID2, varID, missval);
258 }
259 }
260 else if (key == "valid_min")
261 {
262 lvalid_min = true;
263 var.valid_min = parameter_to_double(value);
264 }
265 else if (key == "valid_max")
266 {
267 lvalid_max = true;
268 var.valid_max = parameter_to_double(value);
269 }
270 else if (key == "ok_min_mean_abs")
271 {
272 var.check_min_mean_abs = true;
273 var.ok_min_mean_abs = parameter_to_double(value);
274 }
275 else if (key == "ok_max_mean_abs")
276 {
277 var.check_max_mean_abs = true;
278 var.ok_max_mean_abs = parameter_to_double(value);
279 }
280 else if (key == "datatype" || key == "type")
281 {
282 const auto datatype = cdo_str_to_datatype(parameter_to_word(value.c_str()));
283 if (datatype != -1) vlistDefVarDatatype(vlistID2, varID, datatype);
284 }
285 else
286 {
287 if (Options::cdoVerbose) cdo_print("Attribute %s:%s not supported!", varname, key.c_str());
288 }
289 // clang-format on
290 }
291
292 if (lvalid_min && lvalid_max) var.checkvalid = true;
293 }
294 else
295 {
296 cdo_print("Variable %s not found in CMOR table!", varname);
297 }
298 }
299 }
300
301 void *
CMOR_lite(void * process)302 CMOR_lite(void *process)
303 {
304 bool delvars = false;
305
306 cdo_initialize(process);
307
308 Options::CMOR_Mode = 1;
309 if (Options::CMOR_Mode) cdiDefGlobal("CMOR_MODE", Options::CMOR_Mode);
310
311 cdo_operator_add("cmorlite", 0, 0, "parameter table name");
312
313 const auto operatorID = cdo_operator_id();
314
315 operator_input_arg(cdo_operator_enter(operatorID));
316
317 if (cdo_operator_argc() < 1) cdo_abort("Too few arguments!");
318
319 bool convert_data = false;
320 if (cdo_operator_argc() == 2)
321 {
322 if (cdo_operator_argv(1) == "convert")
323 convert_data = true;
324 else
325 cdo_abort("Unknown parameter: >%s<", cdo_operator_argv(1).c_str());
326 }
327
328 if (cdo_operator_argc() > 2) cdo_abort("Too many arguments!");
329
330 const auto streamID1 = cdo_open_read(0);
331
332 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
333 auto vlistID2 = vlistDuplicate(vlistID1);
334
335 const auto nvars = vlistNvars(vlistID2);
336 std::vector<CmorVar> vars(nvars);
337
338 if (convert_data)
339 for (int varID = 0; varID < nvars; ++varID) vars[varID].convert = true;
340
341 const auto filename = cdo_operator_argv(0).c_str();
342 auto fp = fopen(filename, "r");
343 if (fp == nullptr) cdo_abort("Open failed on: %s\n", filename);
344
345 PMList pmlist;
346 pmlist.read_cmor_table(fp, filename);
347 fclose(fp);
348
349 applyCmorList(pmlist, nvars, vlistID2, vars);
350
351 VarList varList2;
352 varListInit(varList2, vlistID2);
353
354 for (int varID = 0; varID < nvars; ++varID)
355 if (vars[varID].remove)
356 {
357 delvars = true;
358 break;
359 }
360
361 if (delvars)
362 {
363 vlistClearFlag(vlistID1);
364 vlistClearFlag(vlistID2);
365
366 for (int varID = 0; varID < nvars; varID++)
367 {
368 for (int levID = 0; levID < varList2[varID].nlevels; levID++)
369 {
370 vlistDefFlag(vlistID1, varID, levID, true);
371 vlistDefFlag(vlistID2, varID, levID, true);
372 if (vars[varID].remove)
373 {
374 vlistDefFlag(vlistID1, varID, levID, false);
375 vlistDefFlag(vlistID2, varID, levID, false);
376 }
377 }
378 }
379
380 const auto vlistIDx = vlistCreate();
381 cdo_vlist_copy_flag(vlistIDx, vlistID2);
382
383 vlistDestroy(vlistID2);
384
385 vlistID2 = vlistIDx;
386 if (vlistNvars(vlistID2) == 0) cdo_abort("No variable selected!");
387 }
388
389 for (int varID = 0; varID < nvars; ++varID)
390 {
391 auto &var = vars[varID];
392 if (!var.convert) var.changeunits = false;
393 if (var.changeunits)
394 cdo_convert_units(&var.ut_converter, &var.changeunits, (char *) &var.units, (char *) &var.units_old, var.name);
395 }
396
397 const auto taxisID1 = vlistInqTaxis(vlistID1);
398 const auto taxisID2 = taxisDuplicate(taxisID1);
399 vlistDefTaxis(vlistID2, taxisID2);
400
401 /* vlistPrint(vlistID2);*/
402 const auto streamID2 = cdo_open_write(1);
403 cdo_def_vlist(streamID2, vlistID2);
404
405 auto gridsizemax = vlistGridsizeMax(vlistID1);
406 if (vlistNumber(vlistID1) != CDI_REAL) gridsizemax *= 2;
407 Varray<double> array(gridsizemax);
408
409 int tsID = 0;
410 while (true)
411 {
412 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
413 if (nrecs == 0) break;
414
415 cdo_taxis_copy_timestep(taxisID2, taxisID1);
416
417 cdo_def_timestep(streamID2, tsID);
418
419 cmor_check_init(nvars, vars);
420
421 for (int recID = 0; recID < nrecs; recID++)
422 {
423 int varID, levelID;
424 cdo_inq_record(streamID1, &varID, &levelID);
425
426 auto &var = vars[varID];
427 auto varID2 = varID;
428 auto levelID2 = levelID;
429
430 if (delvars)
431 {
432 if (var.remove) continue;
433
434 if (vlistInqFlag(vlistID1, varID, levelID) == true)
435 {
436 varID2 = vlistFindVar(vlistID2, varID);
437 levelID2 = vlistFindLevel(vlistID2, varID, levelID);
438 }
439 }
440
441 cdo_def_record(streamID2, varID2, levelID2);
442
443 size_t nmiss;
444 cdo_read_record(streamID1, array.data(), &nmiss);
445
446 const auto missval = varList2[varID2].missval;
447 auto gridsize = varList2[varID2].gridsize;
448 if (varList2[varID2].nwpv != CDI_REAL) gridsize *= 2;
449
450 if (nmiss && var.changemissval)
451 {
452 for (size_t i = 0; i < gridsize; ++i)
453 {
454 if (DBL_IS_EQUAL(array[i], var.missval_old)) array[i] = missval;
455 }
456 }
457
458 if (var.lfactor)
459 {
460 for (size_t i = 0; i < gridsize; ++i)
461 {
462 if (!DBL_IS_EQUAL(array[i], missval)) array[i] *= var.factor;
463 }
464 }
465
466 #ifdef HAVE_UDUNITS2
467 if (var.changeunits)
468 {
469 int nerr = 0;
470 for (size_t i = 0; i < gridsize; ++i)
471 {
472 if (!DBL_IS_EQUAL(array[i], missval))
473 {
474 array[i] = cv_convert_double((const cv_converter *) var.ut_converter, array[i]);
475 if (ut_get_status() != UT_SUCCESS) nerr++;
476 }
477 }
478 if (nerr)
479 {
480 cdo_warning("Udunits: Error converting units from [%s] to [%s], parameter: %s", var.units_old, var.units,
481 var.name);
482 var.changeunits = false;
483 }
484 }
485 #endif
486
487 cdo_write_record(streamID2, array.data(), nmiss);
488
489 cmor_check_prep(var, gridsize, missval, array.data());
490 }
491
492 cmor_check_eval(vlistID2, nvars, vars);
493
494 tsID++;
495 }
496
497 cdo_stream_close(streamID2);
498 cdo_stream_close(streamID1);
499
500 #ifdef HAVE_UDUNITS2
501 for (int varID = 0; varID < nvars; varID++)
502 if (vars[varID].changeunits) cdo_convert_free(vars[varID].ut_converter);
503
504 cdo_convert_destroy();
505 #endif
506
507 cdo_finish();
508
509 return 0;
510 }
511