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 #include <algorithm>
9
10 #include <cdi.h>
11 #include "cdo_options.h"
12 #include "util_string.h"
13 #include "cdo_vlist.h"
14 #include "compare.h"
15 #include "cdo_output.h"
16 #include "field.h"
17
18 void
vlist_define_timestep_type(int vlistID,int operfunc)19 vlist_define_timestep_type(int vlistID, int operfunc)
20 {
21 int stepType = -1;
22 // clang-format off
23 if (operfunc == FieldFunc_Mean) stepType = TSTEP_AVG;
24 else if (operfunc == FieldFunc_Avg) stepType = TSTEP_AVG;
25 else if (operfunc == FieldFunc_Sum) stepType = TSTEP_SUM;
26 else if (operfunc == FieldFunc_Range) stepType = TSTEP_RANGE;
27 else if (operfunc == FieldFunc_Min) stepType = TSTEP_MIN;
28 else if (operfunc == FieldFunc_Max) stepType = TSTEP_MAX;
29 // clang-format on
30
31 if (stepType != -1)
32 {
33 const auto nvars = vlistNvars(vlistID);
34 for (int varID = 0; varID < nvars; ++varID) vlistDefVarTsteptype(vlistID, varID, stepType);
35 }
36 }
37
38 double
cdo_zaxis_inq_level(int zaxisID,int levelID)39 cdo_zaxis_inq_level(int zaxisID, int levelID)
40 {
41 const auto zaxistype = zaxisInqType(zaxisID);
42 return zaxisInqLevels(zaxisID, nullptr) ? zaxisInqLevel(zaxisID, levelID) : (zaxistype == ZAXIS_SURFACE) ? 0.0 : levelID + 1.0;
43 }
44
45 int
cdo_zaxis_inq_levels(int zaxisID,double * levels)46 cdo_zaxis_inq_levels(int zaxisID, double *levels)
47 {
48 auto size = zaxisInqLevels(zaxisID, nullptr);
49
50 if (levels)
51 {
52 if (size)
53 zaxisInqLevels(zaxisID, levels);
54 else
55 {
56 size = zaxisInqSize(zaxisID);
57 if (size == 1 && zaxisInqType(zaxisID) == ZAXIS_SURFACE)
58 levels[0] = 0.0;
59 else
60 for (int i = 0; i < size; ++i) levels[i] = i + 1.0;
61 }
62 }
63
64 return size;
65 }
66
67 static void
compare_lat_reg2d(size_t ysize,int gridID1,int gridID2)68 compare_lat_reg2d(size_t ysize, int gridID1, int gridID2)
69 {
70 if (ysize > 1)
71 {
72 Varray<double> yvals1(ysize), yvals2(ysize);
73 const auto ny1 = gridInqYvals(gridID1, &yvals1[0]);
74 const auto ny2 = gridInqYvals(gridID2, &yvals2[0]);
75 if (ny1 == 0 || ny2 == 0) return;
76
77 if (IS_EQUAL(yvals1[0], yvals2[ysize - 1]) && IS_EQUAL(yvals1[ysize - 1], yvals2[0]))
78 {
79 if (yvals1[0] > yvals2[0])
80 cdo_abort("Latitude orientation differ! First grid: N->S; second grid: S->N");
81 else
82 cdo_abort("Latitude orientation differ! First grid: S->N; second grid: N->S");
83 }
84 else
85 {
86 for (size_t i = 0; i < ysize; ++i)
87 if (std::fabs(yvals1[i] - yvals2[i]) > 3.e-5)
88 {
89 cdo_warning("Grid latitudes differ!");
90 break;
91 }
92 }
93 }
94 }
95
96 static void
compare_lon_reg2d(size_t xsize,int gridID1,int gridID2)97 compare_lon_reg2d(size_t xsize, int gridID1, int gridID2)
98 {
99 if (xsize > 1)
100 {
101 Varray<double> xvals1(xsize), xvals2(xsize);
102 const auto nx1 = gridInqXvals(gridID1, &xvals1[0]);
103 const auto nx2 = gridInqXvals(gridID2, &xvals2[0]);
104 if (nx1 == 0 || nx2 == 0) return;
105
106 for (size_t i = 0; i < xsize; ++i)
107 if (std::fabs(xvals1[i] - xvals2[i]) > 3.e-5)
108 {
109 cdo_warning("Grid longitudes differ!");
110 break;
111 }
112 }
113 }
114
115 static void
compare_grid_unstructured(int gridID1,int gridID2)116 compare_grid_unstructured(int gridID1, int gridID2)
117 {
118 if (gridInqXvals(gridID1, nullptr) && gridInqXvals(gridID1, nullptr) == gridInqXvals(gridID2, nullptr)
119 && gridInqYvals(gridID1, nullptr) && gridInqYvals(gridID1, nullptr) == gridInqYvals(gridID2, nullptr))
120 {
121 const auto gridsize = gridInqSize(gridID1);
122 Varray<double> xvals1(gridsize), yvals1(gridsize), xvals2(gridsize), yvals2(gridsize);
123 gridInqXvals(gridID1, xvals1.data());
124 gridInqYvals(gridID1, yvals1.data());
125 gridInqXvals(gridID2, xvals2.data());
126 gridInqYvals(gridID2, yvals2.data());
127
128 // size_t inc = gridsize > 10000 ? gridsize / 1000 : 1;
129 size_t inc = 1;
130 for (size_t i = 0; i < gridsize; i += inc)
131 if (std::fabs(xvals1[i] - xvals2[i]) > 2.e-5 || std::fabs(yvals1[i] - yvals2[i]) > 2.e-5)
132 {
133 cdo_warning("Geographic location of some grid points differ!");
134 if (Options::cdoVerbose)
135 printf("cell=%zu x1=%g x2=%g y1=%g y2=%g dx=%g dy=%g\n", i + 1, xvals1[i], xvals2[i], yvals1[i], yvals2[i],
136 xvals1[i] - xvals2[i], yvals1[i] - yvals2[i]);
137 break;
138 }
139 }
140 }
141
142 void
cdo_compare_grids(int gridID1,int gridID2)143 cdo_compare_grids(int gridID1, int gridID2)
144 {
145 if (gridID1 == gridID2) return;
146
147 // compare grids of first variable
148
149 const auto gridType1 = gridInqType(gridID1);
150 const auto gridType2 = gridInqType(gridID2);
151 if (gridType1 == gridType2)
152 {
153 if (gridType1 == GRID_GAUSSIAN || gridType2 == GRID_LONLAT)
154 {
155 const auto xsize = gridInqXsize(gridID1);
156 const auto ysize = gridInqYsize(gridID1);
157
158 if (ysize == gridInqYsize(gridID2))
159 compare_lat_reg2d(ysize, gridID1, gridID2);
160 else
161 cdo_warning("ysize of input grids differ!");
162
163 if (xsize == gridInqXsize(gridID2))
164 compare_lon_reg2d(xsize, gridID1, gridID2);
165 else
166 cdo_warning("xsize of input grids differ!");
167 }
168 else if (gridType1 == GRID_CURVILINEAR || gridType2 == GRID_UNSTRUCTURED)
169 {
170 compare_grid_unstructured(gridID1, gridID2);
171 }
172 }
173 else if (gridInqSize(gridID1) > 1)
174 {
175 cdo_warning("Grids have different types! First grid: %s; second grid: %s", gridNamePtr(gridType1), gridNamePtr(gridType2));
176 }
177 }
178
179 static int
vlistCompareNames(int vlistID1,int varID1,int vlistID2,int varID2)180 vlistCompareNames(int vlistID1, int varID1, int vlistID2, int varID2)
181 {
182 char name1[CDI_MAX_NAME], name2[CDI_MAX_NAME];
183 vlistInqVarName(vlistID1, varID1, name1);
184 vlistInqVarName(vlistID2, varID2, name2);
185 cstr_to_lower_case(name1);
186 cstr_to_lower_case(name2);
187 return strcmp(name1, name2);
188 }
189
190 static int
zaxisCheckLevels(int zaxisID1,int zaxisID2)191 zaxisCheckLevels(int zaxisID1, int zaxisID2)
192 {
193 if (zaxisID1 != zaxisID2)
194 {
195 const auto nlev1 = zaxisInqSize(zaxisID1);
196 const auto nlev2 = zaxisInqSize(zaxisID2);
197 if (nlev1 != nlev2) cdo_abort("Number of levels of the input parameters do not match!");
198
199 Varray<double> lev1(nlev1), lev2(nlev1);
200 cdo_zaxis_inq_levels(zaxisID1, &lev1[0]);
201 cdo_zaxis_inq_levels(zaxisID2, &lev2[0]);
202
203 auto ldiffer = false;
204 for (int i = 0; i < nlev1; ++i)
205 if (IS_NOT_EQUAL(lev1[i], lev2[i]))
206 {
207 ldiffer = true;
208 break;
209 }
210 if (ldiffer)
211 {
212 ldiffer = false;
213 for (int i = 0; i < nlev1; ++i)
214 if (IS_NOT_EQUAL(lev1[i], lev2[nlev1 - 1 - i]))
215 {
216 ldiffer = true;
217 break;
218 }
219
220 if (ldiffer)
221 cdo_warning("Input parameters have different levels!");
222 else
223 cdo_warning("Z-axis orientation differ!");
224
225 return 1;
226 }
227 }
228
229 return 0;
230 }
231
232 static void
vlistCheckNames(int vlistID1,int vlistID2)233 vlistCheckNames(int vlistID1, int vlistID2)
234 {
235 int varID;
236 auto nvars = vlistNvars(vlistID1);
237
238 // std::vector<std::array<char,CDI_MAX_NAME>> names1(nvars); C++14?
239 // std::vector<std::array<char,CDI_MAX_NAME>> names2(nvars);
240 std::vector<std::vector<char>> names1(nvars), names2(nvars);
241 for (varID = 0; varID < nvars; varID++) names1[varID].resize(CDI_MAX_NAME);
242 for (varID = 0; varID < nvars; varID++) names2[varID].resize(CDI_MAX_NAME);
243 for (varID = 0; varID < nvars; varID++) vlistInqVarName(vlistID1, varID, names1[varID].data());
244 for (varID = 0; varID < nvars; varID++) vlistInqVarName(vlistID2, varID, names2[varID].data());
245
246 std::sort(names1.begin(), names1.end());
247 std::sort(names2.begin(), names2.end());
248
249 for (varID = 0; varID < nvars; varID++)
250 if (!cdo_cmpstr(names1[varID].data(), names2[varID].data())) break;
251
252 if (varID == nvars) cdo_print("Use CDO option --sortname to sort the parameter by name (NetCDF only)!");
253 }
254
255 void
vlistPrintMissingVars(int vlistID1,int vlistID2)256 vlistPrintMissingVars(int vlistID1, int vlistID2)
257 {
258 int varID1, varID2;
259 const auto nvars1 = vlistNvars(vlistID1);
260 const auto nvars2 = vlistNvars(vlistID2);
261
262 std::vector<std::vector<char>> names1(nvars1), names2(nvars2);
263 for (varID1 = 0; varID1 < nvars1; varID1++) names1[varID1].resize(CDI_MAX_NAME);
264 for (varID2 = 0; varID2 < nvars2; varID2++) names2[varID2].resize(CDI_MAX_NAME);
265 for (varID1 = 0; varID1 < nvars1; varID1++) vlistInqVarName(vlistID1, varID1, names1[varID1].data());
266 for (varID2 = 0; varID2 < nvars2; varID2++) vlistInqVarName(vlistID2, varID2, names2[varID2].data());
267
268 if (nvars1 > nvars2)
269 {
270 for (varID1 = 0; varID1 < nvars1; varID1++)
271 {
272 for (varID2 = 0; varID2 < nvars2; varID2++)
273 {
274 if (cdo_cmpstr(names1[varID1].data(), names2[varID2].data())) break;
275 }
276 if (varID2 == nvars2) cdo_print("Variable %s not found in second input stream!", names1[varID1].data());
277 }
278 }
279 else
280 {
281 for (varID2 = 0; varID2 < nvars2; varID2++)
282 {
283 for (varID1 = 0; varID1 < nvars1; varID1++)
284 {
285 if (cdo_cmpstr(names1[varID1].data(), names2[varID2].data())) break;
286 }
287 if (varID1 == nvars1) cdo_print("Variable %s not found in first input stream!", names2[varID2].data());
288 }
289 }
290 }
291
292 void
vlist_compare(int vlistID1,int vlistID2,int flag)293 vlist_compare(int vlistID1, int vlistID2, int flag)
294 {
295 auto lchecknames = false;
296
297 const auto nvars = vlistNvars(vlistID1);
298
299 if (nvars != vlistNvars(vlistID2))
300 {
301 vlistPrintMissingVars(vlistID1, vlistID2);
302 cdo_abort("Input streams have different number of variables per timestep!");
303 }
304
305 if (vlistNrecs(vlistID1) != vlistNrecs(vlistID2))
306 cdo_abort("Input streams have different number of %s per timestep!", nvars == 1 ? "layers" : "records");
307
308 for (int varID = 0; varID < nvars; varID++)
309 {
310 if (nvars > 1)
311 {
312 if (flag & CMP_NAME)
313 {
314 if (vlistCompareNames(vlistID1, varID, vlistID2, varID) != 0)
315 {
316 cdo_warning("Input streams have different parameter names!");
317 lchecknames = true;
318 flag -= CMP_NAME;
319 // break;
320 }
321 }
322 }
323
324 if (flag & CMP_GRIDSIZE)
325 {
326 if (gridInqSize(vlistInqVarGrid(vlistID1, varID)) != gridInqSize(vlistInqVarGrid(vlistID2, varID)))
327 {
328 char name[CDI_MAX_NAME];
329 vlistInqVarName(vlistID1, varID, name);
330 cdo_abort("Grid size of the input parameter %s do not match!", name);
331 }
332 }
333
334 if (flag & CMP_NLEVEL)
335 {
336 const auto zaxisID1 = vlistInqVarZaxis(vlistID1, varID);
337 const auto zaxisID2 = vlistInqVarZaxis(vlistID2, varID);
338 if (zaxisCheckLevels(zaxisID1, zaxisID2) != 0) break;
339 }
340 }
341
342 if (flag & CMP_GRID)
343 {
344 const auto gridID1 = vlistInqVarGrid(vlistID1, 0);
345 const auto gridID2 = vlistInqVarGrid(vlistID2, 0);
346 if (gridID1 != gridID2) cdo_compare_grids(gridID1, gridID2);
347 }
348
349 if (lchecknames) vlistCheckNames(vlistID1, vlistID2);
350 }
351
352 int
vlist_compare_x(int vlistID1,int vlistID2,int flag)353 vlist_compare_x(int vlistID1, int vlistID2, int flag)
354 {
355 const auto nvars = vlistNvars(vlistID1);
356 const auto nvars2 = vlistNvars(vlistID2);
357 const auto nlevels2 = zaxisInqSize(vlistInqVarZaxis(vlistID2, 0));
358
359 if (nvars2 != 1) cdo_abort("Internal problem, vlist_compare_x() called with unexpected vlistID2 argument!");
360
361 for (int varID = 0; varID < nvars; varID++)
362 {
363 if (flag & CMP_GRIDSIZE)
364 {
365 if (gridInqSize(vlistInqVarGrid(vlistID1, varID)) != gridInqSize(vlistInqVarGrid(vlistID2, 0)))
366 cdo_abort("Grid size of the input parameters do not match!");
367 }
368
369 if (flag & CMP_NLEVEL)
370 {
371 if ((zaxisInqSize(vlistInqVarZaxis(vlistID1, varID)) != nlevels2) && nlevels2 > 1)
372 cdo_abort("Number of levels of the input parameters do not match!");
373 }
374 }
375
376 if (flag & CMP_GRID)
377 {
378 const auto gridID1 = vlistInqVarGrid(vlistID1, 0);
379 const auto gridID2 = vlistInqVarGrid(vlistID2, 0);
380 if (gridID1 != gridID2) cdo_compare_grids(gridID1, gridID2);
381 }
382
383 return nlevels2;
384 }
385
386 void
vlist_map(int vlistID1,int vlistID2,int flag,int mapflag,std::map<int,int> & mapOfVarIDs)387 vlist_map(int vlistID1, int vlistID2, int flag, int mapflag, std::map<int, int> &mapOfVarIDs)
388 {
389 int varID1, varID2;
390 const auto nvars1 = vlistNvars(vlistID1);
391 const auto nvars2 = vlistNvars(vlistID2);
392
393 std::vector<std::vector<char>> names1(nvars1), names2(nvars2);
394 for (varID1 = 0; varID1 < nvars1; varID1++) names1[varID1].resize(CDI_MAX_NAME);
395 for (varID2 = 0; varID2 < nvars2; varID2++) names2[varID2].resize(CDI_MAX_NAME);
396 for (varID1 = 0; varID1 < nvars1; varID1++) vlistInqVarName(vlistID1, varID1, names1[varID1].data());
397 for (varID2 = 0; varID2 < nvars2; varID2++) vlistInqVarName(vlistID2, varID2, names2[varID2].data());
398
399 if (mapflag == 2)
400 {
401 for (varID2 = 0; varID2 < nvars2; varID2++)
402 {
403 for (varID1 = 0; varID1 < nvars1; varID1++)
404 {
405 if (cdo_cmpstr(names1[varID1].data(), names2[varID2].data())) break;
406 }
407 if (varID1 == nvars1)
408 {
409 cdo_abort("Variable %s not found!", names2[varID2].data());
410 }
411 else
412 {
413 mapOfVarIDs[varID1] = varID2;
414 }
415 }
416 }
417 else
418 {
419 for (varID1 = 0; varID1 < nvars1; varID1++)
420 {
421 for (varID2 = 0; varID2 < nvars2; varID2++)
422 {
423 if (cdo_cmpstr(names1[varID1].data(), names2[varID2].data())) break;
424 }
425 if (varID2 == nvars2)
426 {
427 if (mapflag == 3) continue;
428 cdo_abort("Variable %s not found!", names1[varID1].data());
429 }
430 else
431 {
432 mapOfVarIDs[varID1] = varID2;
433 }
434 }
435 }
436
437 if (mapOfVarIDs.empty()) cdo_abort("No variable found!");
438
439 if (Options::cdoVerbose)
440 for (varID1 = 0; varID1 < nvars1; varID1++)
441 {
442 const auto &it = mapOfVarIDs.find(varID1);
443 if (it != mapOfVarIDs.end())
444 cdo_print("Variable %d:%s mapped to %d:%s", varID1, names1[varID1].data(), it->second, names2[it->second].data());
445 }
446
447 if (mapOfVarIDs.size() > 1)
448 {
449 varID2 = mapOfVarIDs.begin()->second;
450 for (auto it = ++mapOfVarIDs.begin(); it != mapOfVarIDs.end(); ++it)
451 {
452 if (it->second < varID2)
453 cdo_abort("Variable names must be sorted, use CDO option --sortname to sort the parameter by name (NetCDF only)!");
454
455 varID2 = it->second;
456 }
457 }
458
459 for (auto it = mapOfVarIDs.begin(); it != mapOfVarIDs.end(); ++it)
460 {
461 varID1 = it->first;
462 varID2 = it->second;
463
464 if (flag & CMP_GRIDSIZE)
465 {
466 if (gridInqSize(vlistInqVarGrid(vlistID1, varID1)) != gridInqSize(vlistInqVarGrid(vlistID2, varID2)))
467 cdo_abort("Grid size of the input parameters do not match!");
468 }
469
470 if (flag & CMP_NLEVEL)
471 {
472 const auto zaxisID1 = vlistInqVarZaxis(vlistID1, varID1);
473 const auto zaxisID2 = vlistInqVarZaxis(vlistID2, varID2);
474 if (zaxisCheckLevels(zaxisID1, zaxisID2) != 0) break;
475 }
476
477 if (flag & CMP_GRID && varID1 == mapOfVarIDs.begin()->first)
478 {
479 const auto gridID1 = vlistInqVarGrid(vlistID1, varID1);
480 const auto gridID2 = vlistInqVarGrid(vlistID2, varID2);
481 if (gridID1 != gridID2) cdo_compare_grids(gridID1, gridID2);
482 }
483 }
484 }
485
486 bool
vlist_is_szipped(int vlistID)487 vlist_is_szipped(int vlistID)
488 {
489 const auto nvars = vlistNvars(vlistID);
490 for (int varID = 0; varID < nvars; varID++)
491 {
492 const auto comptype = vlistInqVarCompType(vlistID, varID);
493 if (comptype == CDI_COMPRESS_SZIP) return true;
494 }
495
496 return false;
497 }
498
499 int
vlist_inq_nwpv(int vlistID,int varID)500 vlist_inq_nwpv(int vlistID, int varID)
501 {
502 const auto datatype = vlistInqVarDatatype(vlistID, varID);
503 // number of words per value; real:1 complex:2
504 const auto nwpv = (datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64) ? 2 : 1;
505
506 return nwpv;
507 }
508
509 size_t
vlist_check_gridsize(int vlistID)510 vlist_check_gridsize(int vlistID)
511 {
512 auto lerror = false;
513 const auto ngp = gridInqSize(vlistGrid(vlistID, 0));
514
515 // check gridsize
516 const auto ngrids = vlistNgrids(vlistID);
517 for (int index = 0; index < ngrids; ++index)
518 {
519 const auto gridID = vlistGrid(vlistID, index);
520 if (ngp != gridInqSize(gridID))
521 {
522 lerror = true;
523 break;
524 }
525 }
526
527 if (lerror)
528 {
529 cdo_print("This operator requires all variables on the same horizontal grid.");
530 cdo_print("Horizontal grids found:");
531 for (int index = 0; index < ngrids; ++index)
532 {
533 const auto gridID = vlistGrid(vlistID, index);
534 cdo_print(" grid=%d type=%s points=%zu", index + 1, gridNamePtr(gridInqType(gridID)), gridInqSize(gridID));
535 }
536 cdo_abort("The input stream contains variables on different horizontal grids!");
537 }
538
539 return ngp;
540 }
541
542 void
vlist_read_vct(int vlistID,int * rzaxisIDh,int * rnvct,int * rnhlev,int * rnhlevf,int * rnhlevh,Varray<double> & vct)543 vlist_read_vct(int vlistID, int *rzaxisIDh, int *rnvct, int *rnhlev, int *rnhlevf, int *rnhlevh, Varray<double> &vct)
544 {
545 int zaxisIDh = -1;
546 int nhlev = 0, nhlevf = 0, nhlevh = 0;
547 int nvct = 0;
548
549 auto lhavevct = false;
550 auto nzaxis = vlistNzaxis(vlistID);
551 for (int iz = 0; iz < nzaxis; ++iz)
552 {
553 // bool mono_level = false;
554 auto mono_level = true;
555 auto zaxisID = vlistZaxis(vlistID, iz);
556 auto nlevel = zaxisInqSize(zaxisID);
557 auto zaxistype = zaxisInqType(zaxisID);
558
559 if (Options::cdoVerbose)
560 cdo_print("ZAXIS_HYBRID=%d ZAXIS_HYBRID_HALF=%d nlevel=%d mono_level=%d", zaxisInqType(zaxisID) == ZAXIS_HYBRID,
561 zaxisInqType(zaxisID) == ZAXIS_HYBRID_HALF, nlevel, mono_level);
562
563 if ((zaxistype == ZAXIS_HYBRID || zaxistype == ZAXIS_HYBRID_HALF) && nlevel > 1 && !mono_level)
564 {
565 Varray<double> level(nlevel);
566 cdo_zaxis_inq_levels(zaxisID, &level[0]);
567 int l;
568 for (l = 0; l < nlevel; l++)
569 {
570 if ((l + 1) != (int) (level[l] + 0.5)) break;
571 }
572 if (l == nlevel) mono_level = true;
573 }
574
575 if ((zaxistype == ZAXIS_HYBRID || zaxistype == ZAXIS_HYBRID_HALF) && nlevel > 1 && mono_level)
576 {
577 nvct = zaxisInqVctSize(zaxisID);
578 if (nlevel == (nvct / 2 - 1))
579 {
580 if (!lhavevct)
581 {
582 lhavevct = true;
583 zaxisIDh = zaxisID;
584 nhlev = nlevel;
585 nhlevf = nhlev;
586 nhlevh = nhlevf + 1;
587
588 vct.resize(nvct);
589 zaxisInqVct(zaxisID, vct.data());
590 if (Options::cdoVerbose)
591 cdo_print("Detected half-level model definition : nlevel == (nvct/2 - 1) (nlevel: %d, nvct: %d, nhlevf: %d, "
592 "nhlevh: %d) ",
593 nlevel, nvct, nhlevf, nhlevh);
594 }
595 }
596 else if (nlevel == (nvct / 2))
597 {
598 if (!lhavevct)
599 {
600 lhavevct = true;
601 zaxisIDh = zaxisID;
602 nhlev = nlevel;
603 nhlevf = nhlev - 1;
604 nhlevh = nhlev;
605
606 vct.resize(nvct);
607 zaxisInqVct(zaxisID, vct.data());
608 if (Options::cdoVerbose)
609 cdo_print(
610 "Detected full-level model definition : nlevel == (nvct/2) (nlevel: %d, nvct: %d, nhlevf: %d, nhlevh: %d) ",
611 nlevel, nvct, nhlevf, nhlevh);
612 }
613 }
614 else if (nlevel == (nvct - 4 - 1))
615 {
616 if (!lhavevct)
617 {
618 Varray<double> rvct(nvct);
619 zaxisInqVct(zaxisID, rvct.data());
620
621 int voff = 4;
622 if ((int) (rvct[0] + 0.5) == 100000 && rvct[voff] < rvct[voff + 1])
623 {
624 lhavevct = true;
625 zaxisIDh = zaxisID;
626 nhlev = nlevel;
627 nhlevf = nhlev;
628 nhlevh = nhlev + 1;
629
630 int vctsize = 2 * nhlevh;
631 vct.resize(vctsize);
632
633 // calculate VCT for LM
634
635 for (int i = 0; i < vctsize / 2; i++)
636 {
637 if (rvct[voff + i] >= rvct[voff] && rvct[voff + i] <= rvct[3])
638 {
639 vct[i] = rvct[0] * rvct[voff + i];
640 vct[vctsize / 2 + i] = 0;
641 }
642 else
643 {
644 vct[i] = (rvct[0] * rvct[3] * (1 - rvct[voff + i])) / (1 - rvct[3]);
645 vct[vctsize / 2 + i] = (rvct[voff + i] - rvct[3]) / (1 - rvct[3]);
646 }
647 }
648
649 if (Options::cdoVerbose)
650 {
651 for (int i = 0; i < vctsize / 2; i++)
652 fprintf(stdout, "%5d %25.17f %25.17f\n", i, vct[i], vct[vctsize / 2 + i]);
653 }
654 }
655 }
656 }
657 }
658 }
659
660 *rzaxisIDh = zaxisIDh;
661 *rnvct = nvct;
662 *rnhlev = nhlev;
663 *rnhlevf = nhlevf;
664 *rnhlevh = nhlevh;
665 }
666
667 void
vlist_change_hybrid_zaxis(int vlistID1,int vlistID2,int zaxisID1,int zaxisID2)668 vlist_change_hybrid_zaxis(int vlistID1, int vlistID2, int zaxisID1, int zaxisID2)
669 {
670 int nvct0 = 0;
671 Varray<double> vct;
672
673 const auto nzaxis = vlistNzaxis(vlistID1);
674 for (int i = 0; i < nzaxis; ++i)
675 {
676 const auto zaxisID = vlistZaxis(vlistID1, i);
677 const auto nlevel = zaxisInqSize(zaxisID);
678
679 if (zaxisID == zaxisID1 && nlevel > 1)
680 {
681 const auto nvct = zaxisInqVctSize(zaxisID);
682 if (!vct.size())
683 {
684 nvct0 = nvct;
685 vct.resize(nvct);
686 zaxisInqVct(zaxisID, vct.data());
687
688 vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
689 }
690 else
691 {
692 if (nvct0 == nvct && memcmp(vct.data(), zaxisInqVctPtr(zaxisID), nvct * sizeof(double)) == 0)
693 vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
694 }
695 }
696 }
697 }
698
699 int
vlist_get_psvarid(int vlistID,int zaxisID)700 vlist_get_psvarid(int vlistID, int zaxisID)
701 {
702 char psname[CDI_MAX_NAME];
703 int length = CDI_MAX_NAME;
704 cdiInqKeyString(zaxisID, CDI_GLOBAL, CDI_KEY_PSNAME, psname, &length);
705
706 if (psname[0])
707 {
708 char name[CDI_MAX_NAME];
709 const auto nvars = vlistNvars(vlistID);
710 for (int varID = 0; varID < nvars; ++varID)
711 {
712 vlistInqVarName(vlistID, varID, name);
713 if (cdo_cmpstr(name, psname)) return varID;
714 }
715 if (Options::cdoVerbose) cdo_warning("Surface pressure variable not found - %s", psname);
716 }
717
718 return -1;
719 }
720
721 int
vlist_get_first_spectral_grid(int vlistID)722 vlist_get_first_spectral_grid(int vlistID)
723 {
724 // find first spectral grid
725 const auto ngrids = vlistNgrids(vlistID);
726 for (int index = 0; index < ngrids; index++)
727 {
728 const auto gridID = vlistGrid(vlistID, index);
729 if (gridInqType(gridID) == GRID_SPECTRAL) return gridID;
730 }
731
732 return -1;
733 }
734
735 int
vlist_get_first_gaussian_grid(int vlistID)736 vlist_get_first_gaussian_grid(int vlistID)
737 {
738 // find first gaussian grid
739 const auto ngrids = vlistNgrids(vlistID);
740 for (int index = 0; index < ngrids; index++)
741 {
742 const auto gridID = vlistGrid(vlistID, index);
743 if (gridInqType(gridID) == GRID_GAUSSIAN) return gridID;
744 }
745
746 return -1;
747 }
748
749 int
vlist_get_first_fourier_grid(int vlistID)750 vlist_get_first_fourier_grid(int vlistID)
751 {
752 // find first fourier grid
753 const auto ngrids = vlistNgrids(vlistID);
754 for (int index = 0; index < ngrids; index++)
755 {
756 const auto gridID = vlistGrid(vlistID, index);
757 if (gridInqType(gridID) == GRID_FOURIER) return gridID;
758 }
759
760 return -1;
761 }
762
763 void
cdo_check_missval(double missval,const char * varname)764 cdo_check_missval(double missval, const char *varname)
765 {
766 if (DBL_IS_EQUAL(0, missval) || DBL_IS_EQUAL(1, missval))
767 {
768 static bool lwarn = true;
769 if (lwarn)
770 {
771 lwarn = false;
772 cdo_warning("Variable %s has a missing value of %g, this can lead to incorrect results with this operator!", varname,
773 missval);
774 }
775 }
776 }
777