1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Fabian Wachsmann
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 EcaEtccdi eca_etccdi Etccdi conform indices
12 */
13
14 #include <cdi.h>
15
16 #include "cdo_options.h"
17 #include "cdo_vlist.h"
18 #include "datetime.h"
19 #include "process_int.h"
20 #include "param_conversion.h"
21 #include "calendar.h"
22 #include "percentiles_hist.h"
23 #include "percentiles.h"
24 #include "ecautil.h"
25 #include "ecacore.h"
26
27 enum functions_select
28 {
29 func_selle = 1,
30 func_selge = 2
31 };
32
33 static const char TX90P_UNITS[] = "%";
34 static const char TX90P_NAME[] = "tx90pETCCDI";
35 static const char TX90P_LONGNAME[] = "Percentage of Days when Daily Maximum Temperature is Above the 90th Percentile";
36
37 static const char TX10P_UNITS[] = "%";
38 static const char TX10P_NAME[] = "tx10pETCCDI";
39 static const char TX10P_LONGNAME[] = "Percentage of Days when Daily Maximum Temperature is Below the 10th Percentile";
40
41 static const char TN90P_UNITS[] = "%";
42 static const char TN90P_NAME[] = "tn90pETCCDI";
43 static const char TN90P_LONGNAME[] = "Percentage of Days when Daily Minimum Temperature is Above the 90th Percentile";
44
45 static const char TN10P_UNITS[] = "%";
46 static const char TN10P_NAME[] = "tn10pETCCDI";
47 static const char TN10P_LONGNAME[] = "Percentage of Days when Daily Minimum Temperature is Below the 10th Percentile";
48
49 static const char R99P_UNITS[] = "mm";
50 static const char R99P_NAME[] = "r99pETCCDI";
51 static const char R99P_LONGNAME[]
52 = "Annual Total Precipitation when Daily Precipitation Exceeds the 99th Percentile of Wet Day Precipitation";
53
54 static const char R95P_UNITS[] = "mm";
55 static const char R95P_NAME[] = "r95pETCCDI";
56 static const char R95P_LONGNAME[]
57 = "Annual Total Precipitation when Daily Precipitation Exceeds the 95th Percentile of Wet Day Precipitation";
58
59 /* windowDays()
60 Saves for each bootstrap year
61 for each day of year
62 the day of each window day */
63 static void
windowDays(int dayoy,std::vector<int> & wdays,const std::vector<bool> & wdaysRead,int MaxDays,int ndates,int sumboot)64 windowDays(int dayoy, std::vector<int> &wdays, const std::vector<bool> &wdaysRead, int MaxDays, int ndates, int sumboot)
65 {
66 int wdayoy = dayoy;
67 for (int gobackDays = ceil(ndates / 2. - 1); gobackDays != 0; gobackDays--)
68 {
69 wdayoy--;
70 if (wdayoy < 1) wdayoy += (MaxDays - 1) * sumboot;
71 while (wdaysRead[wdayoy] == false)
72 {
73 wdayoy--;
74 if (wdayoy == dayoy) cdo_abort("Too less timesteps!");
75 if (wdayoy < 1) wdayoy += (MaxDays - 1) * sumboot;
76 }
77 }
78 int base = (dayoy - 1) * ndates + 1;
79 wdays[base] = wdayoy;
80 int nndates = 1;
81 while (nndates != ndates)
82 {
83 wdayoy++;
84 if (wdayoy > sumboot * (MaxDays - 1)) wdayoy -= sumboot * (MaxDays - 1);
85 if (wdaysRead[wdayoy] != false)
86 {
87 wdays[base + nndates] = wdayoy;
88 nndates++;
89 }
90 }
91 }
92
93 static void
writeTimesteps(int MaxMonths,int recentYear,FieldVector3D & cei,int frequency,int taxisID4,CdoStreamID streamID4,int * otsID,const std::vector<RecordInfo> & recinfo,int maxrecs,const std::vector<int> & tempdpm,int tempdpy,int func2)94 writeTimesteps(int MaxMonths, int recentYear, FieldVector3D &cei, int frequency, int taxisID4, CdoStreamID streamID4, int *otsID,
95 const std::vector<RecordInfo> &recinfo, int maxrecs, const std::vector<int> &tempdpm, int tempdpy, int func2)
96 {
97 if (frequency == 8)
98 {
99 for (int loopmonth = 2; loopmonth < MaxMonths + 2; loopmonth++)
100 {
101 define_mid_of_time(frequency, taxisID4, recentYear, loopmonth, MaxMonths);
102 cdo_def_timestep(streamID4, *otsID);
103 for (int recID = 0; recID < maxrecs; recID++)
104 {
105 if (*otsID && recinfo[recID].lconst) continue;
106 int varIDo = recinfo[recID].varID;
107 int levelIDo = recinfo[recID].levelID;
108 fieldc_div(cei[loopmonth - 2][0][levelIDo], (double) (tempdpm[loopmonth - 2] / 100.0));
109 cdo_def_record(streamID4, varIDo, levelIDo);
110 cdo_write_record(streamID4, cei[loopmonth - 2][0][levelIDo].vec_d.data(), cei[loopmonth - 2][0][levelIDo].nmiss);
111 }
112 (*otsID)++;
113 }
114 }
115 else
116 {
117 define_mid_of_time(frequency, taxisID4, recentYear, 0, MaxMonths);
118 cdo_def_timestep(streamID4, *otsID);
119 for (int recID = 0; recID < maxrecs; recID++)
120 {
121 if (*otsID && recinfo[recID].lconst) continue;
122 int varIDo = recinfo[recID].varID;
123 int levelIDo = recinfo[recID].levelID;
124 if (func2 == FieldFunc_Avg) fieldc_div(cei[0][0][levelIDo], (double) (tempdpy / 100.0));
125 cdo_def_record(streamID4, varIDo, levelIDo);
126 cdo_write_record(streamID4, cei[0][0][levelIDo].vec_d.data(), cei[0][0][levelIDo].nmiss);
127 }
128 (*otsID)++;
129 }
130 }
131
132 static void
calculateOuterPeriod(Field & field,int MaxMonths,int recentYear,int endOfCalc,FieldVector3D & cei,FieldVector3D & varsPtemp,int frequency,int taxisID4,CdoStreamID streamID4,int * otsID,int vlistID1,const std::vector<RecordInfo> & recinfo,int selection,int func2)133 calculateOuterPeriod(Field &field, int MaxMonths, int recentYear, int endOfCalc, FieldVector3D &cei, FieldVector3D &varsPtemp,
134 int frequency, int taxisID4, CdoStreamID streamID4, int *otsID, int vlistID1,
135 const std::vector<RecordInfo> &recinfo, int selection, int func2)
136 {
137 if (cdo_assert_files_only() == false) cdo_abort("infile1 cannot be a pipe");
138 const auto maxrecs = vlistNrecs(vlistID1);
139 auto cdiStream = streamOpenRead(cdo_get_stream_name(0));
140
141 const auto cdiVlistID = streamInqVlist(cdiStream);
142 const auto cdiTaxisID = vlistInqTaxis(cdiVlistID);
143 std::vector<int> tempdpm(MaxMonths);
144 int tempdpy = 0;
145 for (int i = 0; i < MaxMonths; i++) tempdpm[i] = 0;
146 int year, month, day, tsID = 0, varID, levelID;
147 bool lHasStarted = false;
148
149 if (Options::cdoVerbose) cdo_print("Start to process variables");
150
151 while (true)
152 {
153 const auto nrecs = streamInqTimestep(cdiStream, tsID++);
154 if (nrecs == 0) break;
155
156 auto vdate = taxisInqVdate(cdiTaxisID);
157 cdiDecodeDate(vdate, &year, &month, &day);
158 if (!lHasStarted && year != recentYear)
159 continue;
160 else if (!lHasStarted)
161 lHasStarted = true;
162
163 if (year != recentYear)
164 {
165 writeTimesteps(MaxMonths, recentYear, cei, frequency, taxisID4, streamID4, otsID, recinfo, maxrecs, tempdpm, tempdpy,
166 func2);
167 recentYear = year;
168 tempdpy = 0;
169 tempdpm[0] = 0;
170 field_fill(cei[0][0][0], 0.);
171 if (frequency == 8)
172 for (int loopmonth = 1; loopmonth < MaxMonths; loopmonth++)
173 {
174 tempdpm[loopmonth] = 0;
175 field_fill(cei[loopmonth][0][0], 0.);
176 }
177 }
178 if (year == endOfCalc && func2 == FieldFunc_Avg) break;
179 tempdpy++;
180 auto dayoy = decode_day_of_year(vdate);
181 tempdpm[month - 1]++;
182
183 if (func2 == FieldFunc_Sum) dayoy = 1;
184
185 for (int recID = 0; recID < nrecs; recID++)
186 {
187 streamInqRecord(cdiStream, &varID, &levelID);
188 streamReadRecord(cdiStream, field.vec_d.data(), &field.nmiss);
189
190 Field &pctls = varsPtemp[dayoy][0][levelID];
191 if (selection == func_selle)
192 vfarselle(field, pctls);
193 else if (selection == func_selge)
194 vfarselge(field, pctls);
195
196 auto &array = field.vec_d;
197 if (func2 == FieldFunc_Avg)
198 for (size_t i = 0; i < field.size; i++) array[i] = (DBL_IS_EQUAL(array[i], field.missval)) ? 0.0 : 1.0;
199 else
200 for (size_t i = 0; i < field.size; i++) array[i] = (DBL_IS_EQUAL(array[i], field.missval)) ? 0.0 : array[i];
201 if (frequency == 8)
202 field2_add(cei[(int) ((dayoy - 1) / 31.)][0][levelID], field);
203 else
204 field2_add(cei[0][0][levelID], field);
205 }
206 }
207 if (Options::cdoVerbose) cdo_print("Finished Processing variables");
208 if (year != endOfCalc)
209 writeTimesteps(MaxMonths, year, cei, frequency, taxisID4, streamID4, otsID, recinfo, maxrecs, tempdpm, tempdpy, func2);
210
211 field_fill(cei[0][0][0], 0.);
212
213 if (frequency == 8)
214 for (int loopmonth = 1; loopmonth < MaxMonths; loopmonth++)
215 {
216 tempdpm[loopmonth] = 0;
217 field_fill(cei[loopmonth][0][0], 0.);
218 }
219
220 streamClose(cdiStream);
221 }
222
223 void
etccdi_op(ETCCDI_REQUEST & request)224 etccdi_op(ETCCDI_REQUEST &request)
225 {
226 constexpr int MaxDays = 373;
227 constexpr int MaxMonths = 12;
228 FieldVector2D vars2[MaxDays];
229 HistogramSet hsets[MaxDays];
230
231 const int operatorID = cdo_operator_id();
232 auto selection = cdo_operator_f1(operatorID);
233 auto frequency = cdo_operator_f2(operatorID);
234
235 percentile_set_method("rtype8");
236
237 int FIELD_MEMTYPE = (Options::CDO_Memtype == MemType::Float) ? FIELD_FLT : 0;
238
239 bool wdaysSrc[MaxDays];
240
241 if (request.endboot < request.startboot)
242 {
243 cdo_warning("Your interval end '%d' is before the interval start '%d'. Switched interval years.", request.endboot,
244 request.startboot);
245 request.startboot = request.endboot;
246 request.endboot = request.startboot;
247 }
248 int sumboot = request.endboot - request.startboot + 1;
249 std::vector<bool> wdaysRead((MaxDays - 1) * sumboot + 1);
250 std::vector<int> wdays(request.ndates * (MaxDays - 1) * sumboot + 1);
251 std::vector<int> dpy(sumboot), dpm(MaxMonths * sumboot);
252
253 for (int year = 0; year < sumboot; year++)
254 {
255 dpy[year] = 0;
256 for (int month = 0; month < MaxMonths; month++) dpm[month + year * MaxMonths] = 0;
257 }
258
259 const auto streamID1 = cdo_open_read(0);
260 const auto streamID2 = cdo_open_read(1);
261 const auto streamID3 = cdo_open_read(2);
262
263 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
264 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
265 const auto vlistID3 = cdo_stream_inq_vlist(streamID3);
266 const auto vlistID4 = vlistDuplicate(vlistID1);
267
268 vlist_compare(vlistID1, vlistID2, CMP_ALL);
269 vlist_compare(vlistID1, vlistID3, CMP_ALL);
270
271 const auto taxisID1 = vlistInqTaxis(vlistID1);
272 const auto taxisID2 = vlistInqTaxis(vlistID2);
273 const auto taxisID3 = vlistInqTaxis(vlistID3);
274 // TODO - check that time axes 2 and 3 are equal
275
276 const auto taxisID4 = taxisDuplicate(taxisID1);
277 if (taxisHasBounds(taxisID4)) taxisDeleteBounds(taxisID4);
278 vlistDefTaxis(vlistID4, taxisID4);
279
280 const auto streamID4 = cdo_open_write(3);
281
282 const auto nvars = vlistNvars(vlistID1);
283 bool lOnlyOneVar = true;
284
285 if (nvars == 1)
286 {
287 cdiDefKeyString(vlistID4, 0, CDI_KEY_NAME, request.name);
288 cdiDefKeyString(vlistID4, 0, CDI_KEY_LONGNAME, request.longname);
289 cdiDefKeyString(vlistID4, 0, CDI_KEY_UNITS, request.units);
290 cdiDefAttTxt(vlistID4, 0, "cell_methods", (int) strlen("time: maximum"), "time: maximum");
291 }
292 else
293 {
294 lOnlyOneVar = false;
295 cdo_warning("Your input file has more than one variable. No attributes can be set.");
296 }
297
298 cdo_def_vlist(streamID4, vlistID4);
299
300 const auto maxrecs = vlistNrecs(vlistID1);
301 std::vector<RecordInfo> recinfo(maxrecs);
302
303 const auto gridsizemax = vlistGridsizeMax(vlistID1);
304
305 Field field;
306
307 VarList varList1;
308 varListInit(varList1, vlistID1);
309
310 FieldVector3D vars1(sumboot * (MaxDays - 1) + 1), cei(sumboot * MaxMonths), varsPtemp(MaxDays);
311 // int64_t vdates2[sumboot*(MaxDays-1)+1], vdates1[MaxDays];
312 // int vtimes2[sumboot*(MaxDays-1)+1] , vtimes2[MaxDays];
313 for (int dayoy = 0; dayoy < MaxDays; dayoy++)
314 {
315 fields_from_vlist(vlistID1, vars1[dayoy], FIELD_VEC | FIELD_MEMTYPE);
316 wdaysSrc[dayoy] = false;
317 for (int year = 0; year < sumboot; year++) wdaysRead[dayoy + year * (MaxDays - 1)] = false;
318 }
319
320 for (int dayoy = MaxDays; dayoy < sumboot * (MaxDays - 1) + 1; dayoy++)
321 fields_from_vlist(vlistID1, vars1[dayoy], FIELD_VEC | FIELD_MEMTYPE);
322
323 int tsID = 0;
324
325 while (true)
326 {
327 const auto nrecs = cdo_stream_inq_timestep(streamID2, tsID);
328 if (nrecs == 0) break;
329
330 if (nrecs != cdo_stream_inq_timestep(streamID3, tsID))
331 cdo_abort("Number of records at time step %d of %s and %s differ!", tsID + 1, cdo_get_stream_name(1),
332 cdo_get_stream_name(2));
333
334 auto vdate = taxisInqVdate(taxisID2);
335 // auto vtime = taxisInqVtime(taxisID2);
336
337 if (vdate != taxisInqVdate(taxisID3))
338 cdo_abort("Verification dates at time step %d of %s and %s differ!", tsID + 1, cdo_get_stream_name(1),
339 cdo_get_stream_name(2));
340
341 // if (Options::cdoVerbose) cdo_print("process timestep: %d %d %d", tsID + 1, vdate, vtime);
342
343 auto dayoy = decode_day_of_year(vdate);
344
345 if (request.func2 == FieldFunc_Sum) dayoy = 1;
346
347 if (dayoy < 0 || dayoy >= MaxDays) cdo_abort("Day %d out of range!", dayoy);
348
349 if (!vars2[dayoy].size())
350 {
351 wdaysSrc[dayoy] = true;
352 fields_from_vlist(vlistID2, vars2[dayoy], FIELD_VEC | FIELD_MEMTYPE);
353 fields_from_vlist(vlistID2, varsPtemp[dayoy], FIELD_VEC);
354 hsets[dayoy].create(nvars);
355
356 for (int varID = 0; varID < nvars; varID++)
357 hsets[dayoy].createVarLevels(varID, varList1[varID].nlevels, varList1[varID].gridsize);
358 }
359
360 for (int recID = 0; recID < nrecs; recID++)
361 {
362 int varID, levelID;
363 cdo_inq_record(streamID2, &varID, &levelID);
364 cdo_read_record(streamID2, vars2[dayoy][varID][levelID]);
365 varsPtemp[dayoy][varID][levelID].nmiss = vars2[dayoy][varID][levelID].nmiss;
366 }
367
368 for (int recID = 0; recID < nrecs; recID++)
369 {
370 int varID, levelID;
371 cdo_inq_record(streamID3, &varID, &levelID);
372 field.init(varList1[varID]);
373 cdo_read_record(streamID3, field);
374
375 hsets[dayoy].defVarLevelBounds(varID, levelID, vars2[dayoy][varID][levelID], field);
376 }
377 // fieldsFree(vlistID2, vars2[dayoy]);
378 // fields_from_vlist(vlistID2, vars2[dayoy], FIELD_VEC);
379
380 tsID++;
381 }
382
383 if (Options::cdoVerbose) cdo_print("Defined the boundaries for the histograms");
384 tsID = 0;
385 bool lOnlyRefPeriod = true;
386 int firstYear = 0, lastYear = 0;
387 int year;
388 while (true)
389 {
390 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID++);
391 if (nrecs == 0) break;
392
393 auto vdate = taxisInqVdate(taxisID1);
394 // auto vtime = taxisInqVtime(taxisID1);
395
396 int month, day;
397 cdiDecodeDate(vdate, &year, &month, &day);
398 if (tsID == 1)
399 {
400 if (year > request.startboot)
401 cdo_abort("The interval start year '%d' is before infile start year '%d'.", request.startboot, year);
402 firstYear = year;
403 }
404 lastYear = year;
405
406 if (year >= request.startboot && year <= request.endboot)
407 {
408 auto dayoy = decode_day_of_year(vdate);
409 if (dayoy < 0 || dayoy >= MaxDays) cdo_abort("Day %d out of range!", dayoy);
410 if (wdaysSrc[dayoy] || request.func2 == FieldFunc_Sum)
411 {
412 // Variable independent ?
413 wdaysRead[dayoy + (year - request.startboot) * (MaxDays - 1)] = dayoy + (year - request.startboot) * (MaxDays - 1);
414 dpy[year - request.startboot]++;
415 dpm[(year - request.startboot) * MaxMonths + (int) ((dayoy - 1) / 31.)]++;
416 for (int recID = 0; recID < nrecs; recID++)
417 {
418 int varID, levelID;
419 cdo_inq_record(streamID1, &varID, &levelID);
420 cdo_read_record(streamID1, vars1[dayoy + (year - request.startboot) * (MaxDays - 1)][varID][levelID]);
421 if (tsID == 0)
422 {
423 recinfo[recID].varID = varID;
424 recinfo[recID].levelID = levelID;
425 recinfo[recID].lconst = (varList1[varID].timetype == TIME_CONSTANT);
426 }
427 }
428 }
429 else
430 cdo_warning("Could not find histogram minimum or maximum for day of year: '%d'", dayoy);
431 }
432 else
433 lOnlyRefPeriod = false;
434 }
435 if (Options::cdoVerbose) cdo_print("Read in variables");
436
437 if (year < request.endboot) cdo_abort("The interval end year '%d' is after infile end year '%d'.", request.endboot, year);
438
439 for (year = 0; year < sumboot; year++)
440 {
441 fields_from_vlist(vlistID1, cei[year * MaxMonths], FIELD_VEC);
442 if (frequency == 8)
443 for (int month = 1; month < MaxMonths; month++) fields_from_vlist(vlistID1, cei[year * MaxMonths + month], FIELD_VEC);
444 }
445
446 // printf("Wir beginnen nun mit der Schleife.\n");
447 int bootsyear = 0;
448 int subyear = 0;
449 int otsID = 0;
450
451 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
452 {
453 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
454 {
455 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
456 {
457 windowDays(loopdoy + ytoadd * (MaxDays - 1), wdays, wdaysRead, MaxDays, request.ndates, sumboot);
458 }
459 }
460 }
461 if (Options::cdoVerbose) cdo_print("Calculated window days");
462
463 for (int varID = 0; varID < nvars; varID++)
464 {
465 if (varList1[varID].timetype == TIME_CONSTANT) continue;
466 for (int levelID = 0; levelID < varList1[varID].nlevels; levelID++)
467 {
468 if (request.func2 == FieldFunc_Sum)
469 {
470 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
471 {
472 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
473 {
474 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
475 {
476 auto &source = vars1[loopdoy + ytoadd * (MaxDays - 1)][varID][levelID];
477 auto &hset = hsets[1];
478 hset.addVarLevelValues(varID, levelID, source);
479 }
480 }
481 }
482 }
483 else
484 {
485 #ifdef _OPENMP
486 #pragma omp parallel for shared(sumboot, vars1, request, varID, levelID, hsets, wdays, wdaysRead) schedule(dynamic)
487 #endif
488 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
489 {
490 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
491 {
492 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
493 {
494 for (int ano = 0; ano < request.ndates; ano++)
495 {
496 auto &source = vars1[wdays[ytoadd * request.ndates * (MaxDays - 1) + (loopdoy - 1) * request.ndates
497 + ano + 1]][varID][levelID];
498 auto &hset = hsets[loopdoy];
499 hset.addVarLevelValues(varID, levelID, source);
500 }
501 }
502 }
503 }
504 }
505 }
506 }
507
508 tsID = 0;
509 if (Options::cdoVerbose) cdo_print("Added 30 years to histograms");
510 if (lOnlyOneVar && ((!lOnlyRefPeriod && firstYear != request.startboot) || request.func2 == FieldFunc_Sum))
511 {
512 for (int varID = 0; varID < nvars; varID++)
513 {
514 if (varList1[varID].timetype == TIME_CONSTANT) continue;
515 for (int levelID = 0; levelID < varList1[varID].nlevels; levelID++)
516 {
517 if (request.func2 == FieldFunc_Sum)
518 {
519 auto &pctls = varsPtemp[1][varID][levelID];
520 hsets[1].getVarLevelPercentiles(pctls, varID, levelID, request.pn);
521 }
522 else
523 {
524 #ifdef _OPENMP
525 #pragma omp parallel for shared(request, wdaysSrc, varID, levelID, hsets, varsPtemp) schedule(dynamic)
526 #endif
527 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
528 {
529 if (wdaysSrc[loopdoy])
530 {
531 auto &pctls = varsPtemp[loopdoy][varID][levelID];
532 hsets[loopdoy].getVarLevelPercentiles(pctls, varID, levelID, request.pn);
533 }
534 }
535 }
536 }
537 }
538 field.resize(gridsizemax);
539 calculateOuterPeriod(field, MaxMonths, firstYear, request.startboot, cei, varsPtemp, frequency, taxisID4, streamID4, &otsID,
540 vlistID1, recinfo, selection, request.func2);
541 }
542 else if (!lOnlyRefPeriod && firstYear != request.startboot)
543 cdo_warning("Since you have more than one variable in the input file, only the bootstrapping period can be calculated");
544
545 if (request.func2 == FieldFunc_Avg)
546 {
547 for (bootsyear = request.startboot; bootsyear < request.endboot + 1; bootsyear++)
548 {
549 if (Options::cdoVerbose) cdo_print("Bootsyear: %d", bootsyear);
550 for (int varID = 0; varID < nvars; varID++)
551 {
552 if (varList1[varID].timetype == TIME_CONSTANT) continue;
553 for (int levelID = 0; levelID < varList1[varID].nlevels; levelID++)
554 {
555 #ifdef _OPENMP
556 #pragma omp parallel for shared(sumboot, wdaysRead, request, vars1, varID, levelID, hsets, wdays, cei) schedule(dynamic)
557 #endif
558 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
559 {
560 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
561 {
562 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
563 {
564 for (int ano = 0; ano < request.ndates; ano++)
565 {
566 int recentWday
567 = ytoadd * request.ndates * (MaxDays - 1) + (loopdoy - 1) * request.ndates + ano + 1;
568 if ((int((wdays[recentWday] - 1) / (MaxDays - 1)) + request.startboot) == bootsyear)
569 {
570 auto &source = vars1[wdays[recentWday]][varID][levelID];
571 auto &hset = hsets[loopdoy];
572 hset.subVarLevelValues(varID, levelID, source);
573 }
574 // percyear cannot be smaller than request.startboot
575 if ((int((wdays[recentWday] - 1) / (MaxDays - 1)) + request.startboot) == bootsyear - 1)
576 {
577 auto &source = vars1[wdays[recentWday]][varID][levelID];
578 auto &hset = hsets[loopdoy];
579 hset.addVarLevelValues(varID, levelID, source);
580 }
581 }
582 }
583 }
584 }
585 for (subyear = request.startboot; subyear < request.endboot + 1; subyear++)
586 {
587 if (Options::cdoVerbose) cdo_print("Subyear: %d", subyear);
588 if (subyear != bootsyear)
589 {
590 #ifdef _OPENMP
591 #pragma omp parallel for shared(sumboot, request, vars1, varID, levelID, hsets, wdaysRead, varsPtemp, vars2, cei, subyear, \
592 bootsyear, wdays, frequency) schedule(dynamic)
593 #endif
594 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
595 {
596 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
597 {
598 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
599 {
600 for (int ano = 0; ano < request.ndates; ano++)
601 {
602 int recentWday = ytoadd * request.ndates * (MaxDays - 1)
603 + (loopdoy - 1) * request.ndates + ano + 1;
604 if ((int((wdays[recentWday] - 1) / (MaxDays - 1)) + request.startboot) == subyear)
605 {
606 auto &source = vars1[wdays[recentWday]][varID][levelID];
607 auto &hset = hsets[loopdoy];
608 if (hset.addVarLevelValues(varID, levelID, source) == 1)
609 cdo_print("'%d', '%d", loopdoy, wdays[recentWday]);
610 }
611 }
612 }
613 }
614 // printf("Haben es zum temp array addiert.\n");
615
616 /*** Calculate percentile ***/
617 if (wdaysRead[loopdoy + (bootsyear - request.startboot) * (MaxDays - 1)])
618 {
619 auto &pctls = varsPtemp[loopdoy][varID][levelID];
620 hsets[loopdoy].getVarLevelPercentiles(pctls, varID, levelID, request.pn);
621 /*** Compare data with percentile ***/
622 auto &source = vars1[loopdoy + (bootsyear - request.startboot) * (MaxDays - 1)][varID][levelID];
623 auto &toCompare = vars2[loopdoy][varID][levelID];
624 field_copy(source, toCompare);
625 if (selection == func_selle)
626 vfarselle(toCompare, pctls);
627 else if (selection == func_selge)
628 vfarselge(toCompare, pctls);
629 if (request.func2 == FieldFunc_Avg)
630 {
631 auto &array = toCompare.vec_d;
632 for (size_t i = 0; i < toCompare.size; i++)
633 array[i] = (DBL_IS_EQUAL(array[i], toCompare.missval)) ? 0.0 : 1.0;
634 }
635 else
636 {
637 auto &array = toCompare.vec_d;
638 for (size_t i = 0; i < toCompare.size; i++)
639 array[i] = (DBL_IS_EQUAL(array[i], toCompare.missval)) ? 0.0 : array[i];
640 }
641 // printf("Haben ein Percentil berechnet.\n");
642 // Year sum
643 if (frequency == 8)
644 #ifdef _OPENMP
645 #pragma omp critical
646 #endif
647 field2_add(cei[(bootsyear - request.startboot) * MaxMonths + (int) ((loopdoy - 1) / 31.)][varID]
648 [levelID],
649 toCompare);
650 else
651 #ifdef _OPENMP
652 #pragma omp critical
653 #endif
654 field2_add(cei[(bootsyear - request.startboot) * MaxMonths][varID][levelID], toCompare);
655 }
656 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
657 {
658 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
659 {
660 for (int ano = 0; ano < request.ndates; ano++)
661 {
662 int recentWday = ytoadd * request.ndates * (MaxDays - 1)
663 + (loopdoy - 1) * request.ndates + ano + 1;
664 if ((int((wdays[recentWday] - 1) / (MaxDays - 1)) + request.startboot) == subyear)
665 {
666 auto &source = vars1[wdays[recentWday]][varID][levelID];
667 auto &hset = hsets[loopdoy];
668 if (hset.subVarLevelValues(varID, levelID, source) == 1)
669 cdo_print("'%d', '%d", loopdoy, wdays[recentWday]);
670 }
671 }
672 }
673 }
674 }
675 }
676 }
677 if (frequency == 8)
678 {
679 for (int month = 0; month < MaxMonths; month++)
680 if (cei[(bootsyear - request.startboot) * MaxMonths + month][varID][levelID].vec_d.data())
681 {
682 // Divide vars2 to receive average
683 fieldc_div(cei[(bootsyear - request.startboot) * MaxMonths + month][varID][levelID],
684 (double) ((sumboot - 1) * dpm[(bootsyear - request.startboot) * MaxMonths + month] / 100.0));
685 }
686 }
687 else if (cei[(bootsyear - request.startboot) * MaxMonths][varID][levelID].vec_d.data())
688 {
689 fieldc_div(cei[(bootsyear - request.startboot) * MaxMonths][varID][levelID],
690 (double) ((sumboot - 1) * dpy[bootsyear - request.startboot] / 100.0));
691 }
692 }
693 }
694 if (frequency == 8)
695 {
696 for (int month = 2; month < MaxMonths + 2; month++)
697 {
698 define_mid_of_time(frequency, taxisID4, bootsyear, month, MaxMonths);
699 cdo_def_timestep(streamID4, otsID);
700
701 for (int recID = 0; recID < maxrecs; recID++)
702 {
703 if (otsID && recinfo[recID].lconst) continue;
704 auto varIDo = recinfo[recID].varID;
705 auto levelIDo = recinfo[recID].levelID;
706 cdo_def_record(streamID4, varIDo, levelIDo);
707 cdo_write_record(
708 streamID4, cei[(bootsyear - request.startboot) * MaxMonths + (month - 2)][varIDo][levelIDo].vec_d.data(),
709 cei[(bootsyear - request.startboot) * MaxMonths + (month - 2)][varIDo][levelIDo].nmiss);
710 }
711 otsID++;
712 }
713 }
714 else
715 {
716 define_mid_of_time(frequency, taxisID4, bootsyear, 0, MaxMonths);
717 cdo_def_timestep(streamID4, otsID);
718
719 for (int recID = 0; recID < maxrecs; recID++)
720 {
721 if (otsID && recinfo[recID].lconst) continue;
722 auto varIDo = recinfo[recID].varID;
723 auto levelIDo = recinfo[recID].levelID;
724 cdo_def_record(streamID4, varIDo, levelIDo);
725 cdo_write_record(streamID4, cei[(bootsyear - request.startboot) * MaxMonths][varIDo][levelIDo].vec_d.data(),
726 cei[(bootsyear - request.startboot) * MaxMonths][varIDo][levelIDo].nmiss);
727 }
728 otsID++;
729 }
730 // printf("Haben ein Mittel für Jahr '%d' berechnet.\n", bootsyear);
731 }
732 }
733 if (!lOnlyRefPeriod && lOnlyOneVar && lastYear != request.endboot && request.func2 == FieldFunc_Avg)
734 {
735 field_fill(cei[0][0][0], 0.);
736 if (frequency == 8)
737 for (int loopmonth = 1; loopmonth < MaxMonths; loopmonth++) field_fill(cei[loopmonth][0][0], 0.);
738
739 for (int varID = 0; varID < nvars; varID++)
740 {
741 if (varList1[varID].timetype == TIME_CONSTANT) continue;
742 for (int levelID = 0; levelID < varList1[varID].nlevels; levelID++)
743 {
744 #ifdef _OPENMP
745 #pragma omp parallel for shared(request, wdaysRead, varID, levelID, hsets, varsPtemp) schedule(dynamic)
746 #endif
747 for (int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
748 {
749 for (int ytoadd = 0; ytoadd < sumboot; ytoadd++)
750 {
751 if (wdaysRead[loopdoy + ytoadd * (MaxDays - 1)])
752 {
753 for (int ano = 0; ano < request.ndates; ano++)
754 {
755 int recentWday = ytoadd * request.ndates * (MaxDays - 1) + (loopdoy - 1) * request.ndates + ano + 1;
756 // percyear cannot be smaller than request.startboot
757 if ((int((wdays[recentWday] - 1) / (MaxDays - 1)) + request.startboot) == bootsyear - 1)
758 {
759 auto &source = vars1[wdays[recentWday]][varID][levelID];
760 auto &hset = hsets[loopdoy];
761 hset.addVarLevelValues(varID, levelID, source);
762 }
763 }
764 }
765 }
766 if (wdaysSrc[loopdoy])
767 {
768 auto &pctls = varsPtemp[loopdoy][varID][levelID];
769 hsets[loopdoy].getVarLevelPercentiles(pctls, varID, levelID, request.pn);
770 }
771 }
772 }
773 }
774 field.resize(gridsizemax);
775 field.missval = vars1[1][0][0].missval;
776 field.size = vars1[1][0][0].size;
777 field.grid = vars1[1][0][0].grid;
778 calculateOuterPeriod(field, MaxMonths, request.endboot + 1, lastYear + 1, cei, varsPtemp, frequency, taxisID4, streamID4,
779 &otsID, vlistID1, recinfo, selection, request.func2);
780 }
781
782 cdo_stream_close(streamID4);
783 cdo_stream_close(streamID3);
784 cdo_stream_close(streamID2);
785 cdo_stream_close(streamID1);
786 }
787
788 void *
EcaEtccdi(void * process)789 EcaEtccdi(void *process)
790 {
791 cdo_initialize(process);
792
793 if (cdo_operator_argc() < 3) cdo_abort("Too few arguments!");
794 if (cdo_operator_argc() > 4) cdo_abort("Too many arguments!");
795
796 int TX90P, TX10P, TN90P, TN10P, ALLX, R99P, R95P;
797 if (cdo_operator_argc() == 4 && 'm' == cdo_operator_argv(3)[0])
798 {
799 TX90P = cdo_operator_add("etccdi_tx90p", func_selge, 8, nullptr); /* monthly mode */
800 R99P = cdo_operator_add("etccdi_r99p", func_selge, 8, nullptr); /* monthly mode */
801 R95P = cdo_operator_add("etccdi_r95p", func_selge, 8, nullptr); /* monthly mode */
802 TX10P = cdo_operator_add("etccdi_tx10p", func_selle, 8, nullptr); /* monthly mode */
803 TN90P = cdo_operator_add("etccdi_tn90p", func_selge, 8, nullptr); /* monthly mode */
804 TN10P = cdo_operator_add("etccdi_tn10p", func_selle, 8, nullptr); /* monthly mode */
805 ALLX = cdo_operator_add("etccdi", 0, 8, nullptr); /* monthly mode */
806 }
807 else
808 {
809 TX90P = cdo_operator_add("etccdi_tx90p", func_selge, 31, nullptr);
810 R99P = cdo_operator_add("etccdi_r99p", func_selge, 31, nullptr);
811 R95P = cdo_operator_add("etccdi_r95p", func_selge, 31, nullptr);
812 TX10P = cdo_operator_add("etccdi_tx10p", func_selle, 31, nullptr);
813 TN90P = cdo_operator_add("etccdi_tn90p", func_selge, 31, nullptr);
814 TN10P = cdo_operator_add("etccdi_tn10p", func_selle, 31, nullptr);
815 ALLX = cdo_operator_add("etccdi", 0, 31, nullptr);
816 }
817
818 ETCCDI_REQUEST request;
819
820 request.ndates = parameter_to_int(cdo_operator_argv(0));
821 request.startboot = parameter_to_int(cdo_operator_argv(1));
822
823 const auto operatorID = cdo_operator_id();
824 if (operatorID == TX90P || operatorID == TN90P || operatorID == R95P || operatorID == R99P)
825 {
826 if (operatorID == TX90P || operatorID == TN90P)
827 {
828 if (cdo_operator_argc() < 3)
829 cdo_abort("Operator requires at least 3 parameter values, you provided '%d'!", cdo_operator_argc());
830 request.endboot = parameter_to_int(cdo_operator_argv(2));
831 if (operatorID == TX90P)
832 {
833 request.name = TX90P_NAME;
834 request.longname = TX90P_LONGNAME;
835 request.units = TX90P_UNITS;
836 request.func2 = FieldFunc_Avg;
837 }
838 else if (operatorID == TN90P)
839 {
840 request.name = TN90P_NAME;
841 request.longname = TN90P_LONGNAME;
842 request.units = TN90P_UNITS;
843 request.func2 = FieldFunc_Avg;
844 }
845 request.pn = 90;
846 }
847 else
848 {
849 if (cdo_operator_argc() < 2)
850 cdo_abort("Operator requires at least 2 parameter values, you provided '%d'!", cdo_operator_argc());
851 request.ndates = 1;
852 request.startboot = parameter_to_int(cdo_operator_argv(0));
853 request.endboot = parameter_to_int(cdo_operator_argv(1));
854 if (operatorID == R95P)
855 {
856 request.name = R95P_NAME;
857 request.longname = R95P_LONGNAME;
858 request.units = R95P_UNITS;
859 request.pn = 95;
860 request.func2 = FieldFunc_Sum;
861 }
862 else if (operatorID == R99P)
863 {
864 request.name = R99P_NAME;
865 request.longname = R99P_LONGNAME;
866 request.units = R99P_UNITS;
867 request.pn = 99;
868 request.func2 = FieldFunc_Sum;
869 }
870 }
871 }
872 else if (operatorID == TX10P || operatorID == TN10P)
873 {
874 if (cdo_operator_argc() < 3)
875 cdo_abort("Operator requires at least 3 parameter values, you provided '%d'!", cdo_operator_argc());
876 request.endboot = parameter_to_int(cdo_operator_argv(2));
877 if (operatorID == TX10P)
878 {
879 request.name = TX10P_NAME;
880 request.longname = TX10P_LONGNAME;
881 request.units = TX10P_UNITS;
882 request.func2 = FieldFunc_Avg;
883 }
884 else
885 {
886 request.name = TN10P_NAME;
887 request.longname = TN10P_LONGNAME;
888 request.units = TN10P_UNITS;
889 request.func2 = FieldFunc_Avg;
890 }
891 request.pn = 10;
892 }
893
894 etccdi_op(request);
895 /* else
896 EcaEtccdi(-1, ndates, startboot, endboot); */
897
898 cdo_finish();
899
900 return nullptr;
901 }
902