1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Copyright (C) 2006 Brockmann Consult
5
6 Author: Ralf Quast
7 Uwe Schulzweida
8
9 */
10
11 /*
12 This module contains the following operators:
13
14 Ydrunstat ydrunmin Multi-year daily running minimum
15 Ydrunstat ydrunmax Multi-year daily running maximum
16 Ydrunstat ydrunsum Multi-year daily running sum
17 Ydrunstat ydrunmean Multi-year daily running mean
18 Ydrunstat ydrunavg Multi-year daily running average
19 Ydrunstat ydrunvar Multi-year daily running variance
20 Ydrunstat ydrunvar1 Multi-year daily running variance [Normalize by (n-1)]
21 Ydrunstat ydrunstd Multi-year daily running standard deviation
22 Ydrunstat ydrunstd1 Multi-year daily running standard deviation [Normalize by (n-1)]
23 */
24
25 #include <cdi.h>
26
27 #include "process_int.h"
28 #include "param_conversion.h"
29 #include "calendar.h"
30 #include "datetime.h"
31
32 constexpr int MaxDays = 373;
33
34 class YDAY_STATS
35 {
36 public:
37 int64_t vdate[MaxDays] = { 0 };
38 int vtime[MaxDays] = { 0 };
39 int nsets[MaxDays] = { 0 };
40 FieldVector2D vars1[MaxDays];
41 FieldVector2D vars2[MaxDays];
42 int vlist;
43
YDAY_STATS(int vlistID)44 YDAY_STATS(int vlistID) : vlist(vlistID) {}
45 };
46
47 static void
ydstatUpdate(YDAY_STATS & stats,int64_t vdate,int vtime,const FieldVector2D & vars1,const FieldVector2D & vars2,int nsets,int operfunc)48 ydstatUpdate(YDAY_STATS &stats, int64_t vdate, int vtime, const FieldVector2D &vars1, const FieldVector2D &vars2, int nsets,
49 int operfunc)
50 {
51 const auto lvarstd = (vars2.size() > 0);
52
53 const auto dayoy = decode_day_of_year(vdate);
54 if (dayoy < 0 || dayoy >= MaxDays) cdo_abort("Day %d out of range!", dayoy);
55
56 stats.vdate[dayoy] = vdate;
57 stats.vtime[dayoy] = vtime;
58
59 if (!stats.vars1[dayoy].size())
60 {
61 fields_from_vlist(stats.vlist, stats.vars1[dayoy], FIELD_VEC);
62 if (lvarstd) fields_from_vlist(stats.vlist, stats.vars2[dayoy], FIELD_VEC);
63 }
64
65 const auto nvars = vlistNvars(stats.vlist);
66 for (int varID = 0; varID < nvars; varID++)
67 {
68 if (vlistInqVarTimetype(stats.vlist, varID) == TIME_CONSTANT) continue;
69
70 const auto nlevels = zaxisInqSize(vlistInqVarZaxis(stats.vlist, varID));
71 for (int levelID = 0; levelID < nlevels; levelID++)
72 {
73 if (stats.nsets[dayoy] == 0)
74 {
75 stats.vars1[dayoy][varID][levelID].vec_d = vars1[varID][levelID].vec_d;
76 stats.vars1[dayoy][varID][levelID].nmiss = vars1[varID][levelID].nmiss;
77
78 if (lvarstd)
79 {
80 stats.vars2[dayoy][varID][levelID].vec_d = vars2[varID][levelID].vec_d;
81 stats.vars2[dayoy][varID][levelID].nmiss = vars2[varID][levelID].nmiss;
82 }
83 }
84 else
85 {
86 if (lvarstd)
87 {
88 field2_sum(stats.vars1[dayoy][varID][levelID], vars1[varID][levelID]);
89 field2_sum(stats.vars2[dayoy][varID][levelID], vars2[varID][levelID]);
90 }
91 else
92 {
93 field2_function(stats.vars1[dayoy][varID][levelID], vars1[varID][levelID], operfunc);
94 }
95 }
96 }
97 }
98
99 stats.nsets[dayoy] += nsets;
100 }
101
102 static void
ydstatFinalize(YDAY_STATS & stats,int operfunc)103 ydstatFinalize(YDAY_STATS &stats, int operfunc)
104 {
105 const auto lmean = (operfunc == FieldFunc_Mean || operfunc == FieldFunc_Avg);
106 const auto lstd = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Std1);
107 const auto lvarstd = (lstd || operfunc == FieldFunc_Var || operfunc == FieldFunc_Var1);
108 const int divisor = (operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
109
110 auto fieldc_stdvar_func = lstd ? fieldc_std : fieldc_var;
111
112 for (int dayoy = 0; dayoy < MaxDays; dayoy++)
113 if (stats.nsets[dayoy])
114 {
115 const auto nvars = vlistNvars(stats.vlist);
116 for (int varID = 0; varID < nvars; varID++)
117 {
118 if (vlistInqVarTimetype(stats.vlist, varID) == TIME_CONSTANT) continue;
119
120 const auto nlevels = zaxisInqSize(vlistInqVarZaxis(stats.vlist, varID));
121 for (int levelID = 0; levelID < nlevels; levelID++)
122 {
123 auto nsets = stats.nsets[dayoy];
124 auto &rvars1 = stats.vars1[dayoy][varID][levelID];
125
126 if (lmean)
127 {
128 fieldc_div(rvars1, (double) nsets);
129 }
130 else if (lvarstd)
131 {
132 const auto &rvars2 = stats.vars2[dayoy][varID][levelID];
133 fieldc_stdvar_func(rvars1, rvars2, nsets, divisor);
134 }
135 }
136 }
137 }
138 }
139
140 static void
addOperators(void)141 addOperators(void)
142 {
143 // clang-format off
144 cdo_operator_add("ydrunmin", FieldFunc_Min, 0, nullptr);
145 cdo_operator_add("ydrunmax", FieldFunc_Max, 0, nullptr);
146 cdo_operator_add("ydrunsum", FieldFunc_Sum, 0, nullptr);
147 cdo_operator_add("ydrunmean", FieldFunc_Mean, 0, nullptr);
148 cdo_operator_add("ydrunavg", FieldFunc_Avg, 0, nullptr);
149 cdo_operator_add("ydrunvar", FieldFunc_Var, 0, nullptr);
150 cdo_operator_add("ydrunvar1", FieldFunc_Var1, 0, nullptr);
151 cdo_operator_add("ydrunstd", FieldFunc_Std, 0, nullptr);
152 cdo_operator_add("ydrunstd1", FieldFunc_Std1, 0, nullptr);
153 // clang-format on
154 }
155
156 void *
Ydrunstat(void * process)157 Ydrunstat(void *process)
158 {
159 cdo_initialize(process);
160
161 addOperators();
162
163 const auto operatorID = cdo_operator_id();
164 const auto operfunc = cdo_operator_f1(operatorID);
165
166 operator_input_arg("number of timesteps");
167 const auto nparams = cdo_operator_argc();
168 const auto ndates = parameter_to_int(cdo_operator_argv(0));
169 const char readMethod = (nparams == 2) ? cdo_operator_argv(1)[0] : '0';
170
171 const auto lvarstd
172 = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Var || operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
173
174 const auto streamID1 = cdo_open_read(0);
175
176 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
177 const auto vlistID2 = vlistDuplicate(vlistID1);
178
179 const auto taxisID1 = vlistInqTaxis(vlistID1);
180 const auto taxisID2 = taxisDuplicate(taxisID1);
181 if (taxisHasBounds(taxisID2)) taxisDeleteBounds(taxisID2);
182 vlistDefTaxis(vlistID2, taxisID2);
183
184 const auto dpy = calendar_dpy(taxisInqCalendar(taxisID1));
185
186 const auto streamID2 = cdo_open_write(1);
187 cdo_def_vlist(streamID2, vlistID2);
188
189 const auto maxrecs = vlistNrecs(vlistID1);
190 std::vector<RecordInfo> recList(maxrecs);
191
192 std::vector<CdoDateTime> datetime(ndates + 1);
193
194 YDAY_STATS stats(vlistID1);
195 FieldVector3D vars1(ndates + 1), vars2(ndates + 1);
196
197 for (int its = 0; its < ndates; its++)
198 {
199 fields_from_vlist(vlistID1, vars1[its], FIELD_VEC);
200 if (lvarstd) fields_from_vlist(vlistID1, vars2[its], FIELD_VEC);
201 }
202
203 int tsID = 0;
204 int startYear = 0, endYear = 0, mon = 0, day = 0;
205
206 for (tsID = 0; tsID < ndates; tsID++)
207 {
208 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
209 if (nrecs == 0) cdo_abort("File has less then %d timesteps!", ndates);
210
211 datetime[tsID].date = taxisInqVdate(taxisID1);
212 datetime[tsID].time = taxisInqVtime(taxisID1);
213
214 if (tsID == 0 && readMethod == 'c') cdiDecodeDate(datetime[tsID].date, &startYear, &mon, &day);
215
216 for (int recID = 0; recID < nrecs; recID++)
217 {
218 int varID, levelID;
219 cdo_inq_record(streamID1, &varID, &levelID);
220
221 if (tsID == 0)
222 {
223 recList[recID].varID = varID;
224 recList[recID].levelID = levelID;
225 recList[recID].lconst = (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT);
226 }
227
228 auto &rvars1 = vars1[tsID][varID][levelID];
229
230 cdo_read_record(streamID1, rvars1);
231
232 if (lvarstd)
233 {
234 field2_moq(vars2[tsID][varID][levelID], rvars1);
235 for (int inp = 0; inp < tsID; inp++) field2_sumsumq(vars1[inp][varID][levelID], vars2[inp][varID][levelID], rvars1);
236 }
237 else
238 {
239 for (int inp = 0; inp < tsID; inp++) field2_function(vars1[inp][varID][levelID], rvars1, operfunc);
240 }
241 }
242 }
243
244 while (true)
245 {
246 datetime_avg(dpy, ndates, datetime.data());
247
248 const auto vdate = datetime[ndates].date;
249 const auto vtime = datetime[ndates].time;
250
251 ydstatUpdate(stats, vdate, vtime, vars1[0], vars2[0], ndates, operfunc);
252
253 datetime[ndates] = datetime[0];
254 vars1[ndates] = vars1[0];
255 if (lvarstd) vars2[ndates] = vars2[0];
256
257 for (int inp = 0; inp < ndates; inp++)
258 {
259 datetime[inp] = datetime[inp + 1];
260 vars1[inp] = vars1[inp + 1];
261 if (lvarstd) vars2[inp] = vars2[inp + 1];
262 }
263
264 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
265 if (nrecs == 0) break;
266
267 datetime[ndates - 1].date = taxisInqVdate(taxisID1);
268 datetime[ndates - 1].time = taxisInqVtime(taxisID1);
269
270 for (int recID = 0; recID < nrecs; recID++)
271 {
272 int varID, levelID;
273 cdo_inq_record(streamID1, &varID, &levelID);
274
275 auto &rvars1 = vars1[ndates - 1][varID][levelID];
276
277 cdo_read_record(streamID1, rvars1);
278
279 if (lvarstd)
280 {
281 field2_moq(vars2[ndates - 1][varID][levelID], rvars1);
282 for (int inp = 0; inp < ndates - 1; inp++)
283 field2_sumsumq(vars1[inp][varID][levelID], vars2[inp][varID][levelID], rvars1);
284 }
285 else
286 {
287 for (int inp = 0; inp < ndates - 1; inp++) field2_function(vars1[inp][varID][levelID], rvars1, operfunc);
288 }
289 }
290
291 tsID++;
292 }
293
294 if (readMethod == 'c' && cdo_assert_files_only())
295 {
296 cdiDecodeDate(datetime[ndates - 1].date, &endYear, &mon, &day);
297 const auto cdiStream = streamOpenRead(cdo_get_stream_name(0));
298 const auto cdiVlistID = streamInqVlist(cdiStream);
299 const auto cdiTaxisID = vlistInqTaxis(cdiVlistID);
300 int missTimes = 0;
301 for (missTimes = 0; missTimes < ndates - 1; missTimes++)
302 {
303 const auto nrecs = streamInqTimestep(cdiStream, missTimes);
304 if (nrecs == 0) break;
305
306 int year;
307 cdiDecodeDate(taxisInqVdate(cdiTaxisID), &year, &mon, &day);
308 datetime[ndates - 1].date = cdiEncodeDate(endYear + 1, mon, day);
309 datetime[ndates - 1].time = taxisInqVtime(cdiTaxisID);
310
311 for (int recID = 0; recID < nrecs; recID++)
312 {
313 int varID, levelID;
314 streamInqRecord(cdiStream, &varID, &levelID);
315
316 auto &rvars1 = vars1[ndates - 1][varID][levelID];
317
318 streamReadRecord(cdiStream, rvars1.vec_d.data(), &rvars1.nmiss);
319
320 if (lvarstd)
321 {
322 field2_moq(vars2[ndates - 1][varID][levelID], rvars1);
323 for (int inp = 0; inp < ndates - 1; inp++)
324 field2_sumsumq(vars1[inp][varID][levelID], vars2[inp][varID][levelID], rvars1);
325 }
326 else
327 {
328 for (int inp = 0; inp < ndates - 1; inp++) field2_function(vars1[inp][varID][levelID], rvars1, operfunc);
329 }
330 }
331 datetime_avg(dpy, ndates, datetime.data());
332 cdiDecodeDate(datetime[ndates].date, &year, &mon, &day);
333 const auto vdate = (year > endYear) ? cdiEncodeDate(startYear, mon, day) : datetime[ndates].date;
334 const auto vtime = datetime[ndates].time;
335
336 ydstatUpdate(stats, vdate, vtime, vars1[0], vars2[0], ndates, operfunc);
337
338 datetime[ndates] = datetime[0];
339 vars1[ndates] = vars1[0];
340 if (lvarstd) vars2[ndates] = vars2[0];
341
342 for (int inp = 0; inp < ndates; inp++)
343 {
344 datetime[inp] = datetime[inp + 1];
345 vars1[inp] = vars1[inp + 1];
346 if (lvarstd) vars2[inp] = vars2[inp + 1];
347 }
348 }
349
350 if (missTimes != ndates - 1) cdo_abort("Addding the missing values when using the 'readMethod' method was not possible");
351
352 streamClose(cdiStream);
353 }
354 else if (readMethod == 'c')
355 cdo_warning("Operators cannot be piped in circular mode");
356
357 ydstatFinalize(stats, operfunc);
358
359 int otsID = 0;
360
361 for (int dayoy = 0; dayoy < MaxDays; dayoy++)
362 if (stats.nsets[dayoy])
363 {
364 taxisDefVdate(taxisID2, stats.vdate[dayoy]);
365 taxisDefVtime(taxisID2, stats.vtime[dayoy]);
366 cdo_def_timestep(streamID2, otsID);
367
368 for (int recID = 0; recID < maxrecs; recID++)
369 {
370 if (otsID && recList[recID].lconst) continue;
371
372 const auto varID = recList[recID].varID;
373 const auto levelID = recList[recID].levelID;
374 auto &rvars1 = stats.vars1[dayoy][varID][levelID];
375
376 cdo_def_record(streamID2, varID, levelID);
377 cdo_write_record(streamID2, rvars1);
378 }
379
380 otsID++;
381 }
382
383 cdo_stream_close(streamID2);
384 cdo_stream_close(streamID1);
385
386 cdo_finish();
387
388 return nullptr;
389 }
390