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