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 Comp eq Equal
12 Comp ne Not equal
13 Comp le Less equal
14 Comp lt Less than
15 Comp ge Greater equal
16 Comp gt Greater than
17 */
18
19 #include <cdi.h>
20
21 #include <utility>
22
23 #include "cdo_options.h"
24 #include "process_int.h"
25 #include "cdo_vlist.h"
26 #include "cdo_cdi_wrapper.h"
27 #include "cdo_fill.h"
28
29 static
30 auto func_comp = [](auto hasMissvals, auto n, auto mv1, auto mv2, const auto &vIn1, const auto &vIn2, auto &vOut, auto binray_operator)
__anon5497b5f20102(auto hasMissvals, auto n, auto mv1, auto mv2, const auto &vIn1, const auto &vIn2, auto &vOut, auto binray_operator) 31 {
32 if (hasMissvals)
33 {
34 if (std::isnan(mv1) || std::isnan(mv2))
35 for (size_t i = 0; i < n; ++i) vOut[i] = (dbl_is_equal(vIn1[i], mv1) || dbl_is_equal(vIn2[i], mv2)) ? mv1 : binray_operator(vIn1[i], vIn2[i]);
36 else
37 for (size_t i = 0; i < n; ++i) vOut[i] = (is_equal(vIn1[i], mv1) || is_equal(vIn2[i], mv2)) ? mv1 : binray_operator(vIn1[i], vIn2[i]);
38 }
39 else
40 {
41 for (size_t i = 0; i < n; ++i) vOut[i] = binray_operator(vIn1[i], vIn2[i]);
42 }
43 };
44
45 void *
Comp(void * process)46 Comp(void *process)
47 {
48 enum
49 {
50 FILL_NONE,
51 FILL_TS,
52 FILL_REC
53 };
54 int filltype = FILL_NONE;
55 Varray2D<double> vardata;
56
57 cdo_initialize(process);
58
59 const auto EQ = cdo_operator_add("eq", 0, 0, nullptr);
60 const auto NE = cdo_operator_add("ne", 0, 0, nullptr);
61 const auto LE = cdo_operator_add("le", 0, 0, nullptr);
62 const auto LT = cdo_operator_add("lt", 0, 0, nullptr);
63 const auto GE = cdo_operator_add("ge", 0, 0, nullptr);
64 const auto GT = cdo_operator_add("gt", 0, 0, nullptr);
65
66 const auto operatorID = cdo_operator_id();
67
68 operator_check_argc(0);
69
70 auto streamID1 = cdo_open_read(0);
71 auto streamID2 = cdo_open_read(1);
72
73 auto vlistID1 = cdo_stream_inq_vlist(streamID1);
74 auto vlistID2 = cdo_stream_inq_vlist(streamID2);
75
76 auto taxisID1 = vlistInqTaxis(vlistID1);
77 auto taxisID2 = vlistInqTaxis(vlistID2);
78
79 auto ntsteps1 = vlistNtsteps(vlistID1);
80 auto ntsteps2 = vlistNtsteps(vlistID2);
81 if (ntsteps1 == 0) ntsteps1 = 1;
82 if (ntsteps2 == 0) ntsteps2 = 1;
83
84 auto fillstream1 = false;
85
86 if (vlistNrecs(vlistID1) != 1 && vlistNrecs(vlistID2) == 1)
87 {
88 filltype = FILL_REC;
89 cdo_print("Filling up stream2 >%s< by copying the first record.", cdo_get_stream_name(1));
90 if (ntsteps2 != 1) cdo_abort("stream2 has more than 1 timestep!");
91 }
92 else if (vlistNrecs(vlistID1) == 1 && vlistNrecs(vlistID2) != 1)
93 {
94 filltype = FILL_REC;
95 cdo_print("Filling up stream1 >%s< by copying the first record.", cdo_get_stream_name(0));
96 if (ntsteps1 != 1) cdo_abort("stream1 has more than 1 timestep!");
97 fillstream1 = true;
98 std::swap(streamID1, streamID2);
99 std::swap(vlistID1, vlistID2);
100 std::swap(taxisID1, taxisID2);
101 }
102
103 if (filltype == FILL_NONE) vlist_compare(vlistID1, vlistID2, CMP_ALL);
104
105 const auto gridsizemax = vlistGridsizeMax(vlistID1);
106 Varray<double> vaIn1(gridsizemax), vaIn2(gridsizemax), vaOut(gridsizemax);
107
108 double *arrayx1 = vaIn1.data();
109 double *arrayx2 = vaIn2.data();
110
111 if (Options::cdoVerbose) cdo_print("Number of timesteps: file1 %d, file2 %d", ntsteps1, ntsteps2);
112
113 if (filltype == FILL_NONE)
114 {
115 if (ntsteps1 != 1 && ntsteps2 == 1)
116 {
117 filltype = FILL_TS;
118 cdo_print("Filling up stream2 >%s< by copying the first timestep.", cdo_get_stream_name(1));
119 }
120 else if (ntsteps1 == 1 && ntsteps2 != 1)
121 {
122 filltype = FILL_TS;
123 cdo_print("Filling up stream1 >%s< by copying the first timestep.", cdo_get_stream_name(0));
124 fillstream1 = true;
125 std::swap(streamID1, streamID2);
126 std::swap(vlistID1, vlistID2);
127 std::swap(taxisID1, taxisID2);
128 }
129
130 if (filltype == FILL_TS) cdo_fill_ts(vlistID2, vardata);
131 }
132
133 if (fillstream1)
134 {
135 arrayx1 = vaIn2.data();
136 arrayx2 = vaIn1.data();
137 }
138
139 VarList varList1, varList2;
140 varListInit(varList1, vlistID1);
141 varListInit(varList2, vlistID2);
142
143 const auto vlistID3 = vlistDuplicate(vlistID1);
144
145 const auto taxisID3 = taxisDuplicate(taxisID1);
146 vlistDefTaxis(vlistID3, taxisID3);
147
148 const auto streamID3 = cdo_open_write(2);
149 cdo_def_vlist(streamID3, vlistID3);
150
151 int tsID = 0;
152 while (true)
153 {
154 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
155 if (nrecs == 0) break;
156
157 if (tsID == 0 || filltype == FILL_NONE)
158 {
159 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
160 if (nrecs2 == 0) cdo_abort("Input streams have different number of timesteps!");
161 }
162
163 cdo_taxis_copy_timestep(taxisID3, taxisID1);
164 cdo_def_timestep(streamID3, tsID);
165
166 for (int recID = 0; recID < nrecs; recID++)
167 {
168 size_t nmiss1 = 0, nmiss2 = 0;
169 int varID, levelID;
170 cdo_inq_record(streamID1, &varID, &levelID);
171 cdo_read_record(streamID1, arrayx1, &nmiss1);
172
173 if (tsID == 0 || filltype == FILL_NONE)
174 {
175 if (recID == 0 || filltype != FILL_REC)
176 {
177 cdo_inq_record(streamID2, &varID, &levelID);
178 cdo_read_record(streamID2, arrayx2, &nmiss2);
179 }
180
181 if (filltype == FILL_TS)
182 {
183 const auto offset = varList2[varID].gridsize * levelID;
184 array_copy(varList2[varID].gridsize, arrayx2, &vardata[varID][offset]);
185 }
186 }
187 else if (filltype == FILL_TS)
188 {
189 const auto offset = varList2[varID].gridsize * levelID;
190 array_copy(varList2[varID].gridsize, &vardata[varID][offset], arrayx2);
191 }
192
193 const auto datatype1 = varList1[varID].datatype;
194 const auto gridsize1 = varList1[varID].gridsize;
195 auto missval1 = varList1[varID].missval;
196
197 const auto xvarID = (filltype == FILL_REC) ? 0 : varID;
198 const auto datatype2 = varList2[xvarID].datatype;
199 const auto gridsize2 = varList2[xvarID].gridsize;
200 auto missval2 = varList2[xvarID].missval;
201
202 if (gridsize1 != gridsize2)
203 cdo_abort("Streams have different gridsize (gridsize1 = %zu; gridsize2 = %zu)!", gridsize1, gridsize2);
204
205 const auto ngp = gridsize1;
206
207 if (datatype1 != datatype2)
208 {
209 if (datatype1 == CDI_DATATYPE_FLT32 && datatype2 == CDI_DATATYPE_FLT64)
210 {
211 missval2 = (float) missval2;
212 for (size_t i = 0; i < ngp; i++) vaIn2[i] = (float) vaIn2[i];
213 }
214 else if (datatype1 == CDI_DATATYPE_FLT64 && datatype2 == CDI_DATATYPE_FLT32)
215 {
216 missval1 = (float) missval1;
217 for (size_t i = 0; i < ngp; i++) vaIn1[i] = (float) vaIn1[i];
218 }
219 }
220
221 if (nmiss1 > 0) cdo_check_missval(missval1, varList1[varID].name);
222 // if (nmiss2 > 0) cdo_check_missval(missval2, varList2[varID].name);
223
224 const auto hasMissvals = (nmiss1 > 0 || nmiss2 > 0);
225 // clang-format off
226 if (operatorID == EQ) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_EQ);
227 else if (operatorID == NE) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_NE);
228 else if (operatorID == LE) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_LE);
229 else if (operatorID == LT) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_LT);
230 else if (operatorID == GE) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_GE);
231 else if (operatorID == GT) func_comp(hasMissvals, ngp, missval1, missval2, vaIn1, vaIn2, vaOut, binary_op_GT);
232 else cdo_abort("Operator not implemented!");
233 // clang-format on
234
235 const auto nmissOut = varray_num_mv(ngp, vaOut, missval1);
236 cdo_def_record(streamID3, varID, levelID);
237 cdo_write_record(streamID3, vaOut.data(), nmissOut);
238 }
239
240 tsID++;
241 }
242
243 cdo_stream_close(streamID3);
244 cdo_stream_close(streamID2);
245 cdo_stream_close(streamID1);
246
247 cdo_finish();
248
249 return nullptr;
250 }
251