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 Timpctl timpctl Time percentiles
15 Hourpctl hourpctl Hourly percentiles
16 Daypctl daypctl Daily percentiles
17 Monpctl monpctl Monthly percentiles
18 Yearpctl yearpctl Yearly percentiles
19 */
20
21 #include <cdi.h>
22
23 #include "util_date.h"
24 #include "process_int.h"
25 #include "cdo_vlist.h"
26 #include "param_conversion.h"
27 #include "percentiles_hist.h"
28 #include "datetime.h"
29
30 static void
timpctl(int operatorID)31 timpctl(int operatorID)
32 {
33 const auto timestat_date = TimeStat::MEAN;
34 char indate1[DATE_LEN + 1], indate2[DATE_LEN + 1];
35
36 operator_input_arg("percentile number");
37 const auto pn = parameter_to_double(cdo_operator_argv(0));
38
39 const int cmplen = DATE_LEN - cdo_operator_f2(operatorID);
40
41 const auto streamID1 = cdo_open_read(0);
42 const auto streamID2 = cdo_open_read(1);
43 const auto streamID3 = cdo_open_read(2);
44
45 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
46 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
47 const auto vlistID3 = cdo_stream_inq_vlist(streamID3);
48 const auto vlistID4 = vlistDuplicate(vlistID1);
49
50 vlist_compare(vlistID1, vlistID2, CMP_ALL);
51 vlist_compare(vlistID1, vlistID3, CMP_ALL);
52
53 if (cdo_operator_f2(operatorID) == 16) vlistDefNtsteps(vlistID4, 1);
54
55 const auto taxisID1 = vlistInqTaxis(vlistID1);
56 const auto taxisID2 = vlistInqTaxis(vlistID2);
57 const auto taxisID3 = vlistInqTaxis(vlistID3);
58 // TODO - check that time axes 2 and 3 are equal
59
60 const auto taxisID4 = taxisDuplicate(taxisID1);
61 taxisWithBounds(taxisID4);
62 vlistDefTaxis(vlistID4, taxisID4);
63
64 const auto streamID4 = cdo_open_write(3);
65 cdo_def_vlist(streamID4, vlistID4);
66
67 const auto ntsteps = vlistNtsteps(vlistID1);
68 const auto nvars = vlistNvars(vlistID1);
69
70 const auto maxrecs = vlistNrecs(vlistID1);
71 std::vector<RecordInfo> recList(maxrecs);
72
73 FieldVector constFields(maxrecs);
74
75 DateTimeList dtlist;
76 dtlist.set_stat(timestat_date);
77 dtlist.set_calendar(taxisInqCalendar(taxisID1));
78
79 Field field1, field2;
80
81 VarList varList1;
82 varListInit(varList1, vlistID1);
83
84 HistogramSet hset(nvars, ntsteps);
85
86 for (int varID = 0; varID < nvars; varID++) hset.createVarLevels(varID, varList1[varID].nlevels, varList1[varID].gridsize);
87
88 int tsID = 0;
89 int otsID = 0;
90 while (true)
91 {
92 auto nrecs = cdo_stream_inq_timestep(streamID2, otsID);
93 if (nrecs != cdo_stream_inq_timestep(streamID3, otsID))
94 cdo_abort("Number of records at time step %d of %s and %s differ!", otsID + 1, cdo_get_stream_name(1),
95 cdo_get_stream_name(2));
96
97 const auto vdate2 = taxisInqVdate(taxisID2);
98 const auto vtime2 = taxisInqVtime(taxisID2);
99 const auto vdate3 = taxisInqVdate(taxisID3);
100 const auto vtime3 = taxisInqVtime(taxisID3);
101 if (vdate2 != vdate3 || vtime2 != vtime3)
102 cdo_abort("Verification dates at time step %d of %s and %s differ!", otsID + 1, cdo_get_stream_name(1),
103 cdo_get_stream_name(2));
104
105 for (int recID = 0; recID < nrecs; recID++)
106 {
107 int varID, levelID;
108 cdo_inq_record(streamID2, &varID, &levelID);
109 field1.init(varList1[varID]);
110 cdo_read_record(streamID2, field1);
111
112 cdo_inq_record(streamID3, &varID, &levelID);
113 field2.init(varList1[varID]);
114 cdo_read_record(streamID3, field2);
115
116 hset.defVarLevelBounds(varID, levelID, field1, field2);
117 }
118
119 int nsets = 0;
120 while (nrecs && (nrecs = cdo_stream_inq_timestep(streamID1, tsID)))
121 {
122 dtlist.taxis_inq_timestep(taxisID1, nsets);
123 const auto vdate1 = dtlist.get_vdate(nsets);
124 const auto vtime1 = dtlist.get_vtime(nsets);
125
126 if (nsets == 0) SET_DATE(indate2, vdate1, vtime1);
127 SET_DATE(indate1, vdate1, vtime1);
128
129 if (DATE_IS_NEQ(indate1, indate2, cmplen)) break;
130
131 for (int recID = 0; recID < nrecs; recID++)
132 {
133 int varID, levelID;
134 cdo_inq_record(streamID1, &varID, &levelID);
135
136 if (tsID == 0)
137 {
138 recList[recID].varID = varID;
139 recList[recID].levelID = levelID;
140 recList[recID].lconst = (varList1[varID].timetype == TIME_CONSTANT);
141 }
142
143 if (tsID == 0 && recList[recID].lconst)
144 {
145 constFields[recID].init(varList1[varID]);
146 cdo_read_record(streamID1, constFields[recID]);
147 }
148 else
149 {
150 field1.init(varList1[varID]);
151 cdo_read_record(streamID1, field1);
152
153 hset.addVarLevelValues(varID, levelID, field1);
154 }
155 }
156
157 nsets++;
158 tsID++;
159 }
160
161 if (nrecs == 0 && nsets == 0) break;
162
163 dtlist.stat_taxis_def_timestep(taxisID4, nsets);
164 cdo_def_timestep(streamID4, otsID);
165
166 for (int recID = 0; recID < maxrecs; recID++)
167 {
168 if (otsID && recList[recID].lconst) continue;
169
170 const auto varID = recList[recID].varID;
171 const auto levelID = recList[recID].levelID;
172 cdo_def_record(streamID4, varID, levelID);
173
174 if (recList[recID].lconst)
175 {
176 cdo_write_record(streamID4, constFields[recID]);
177 }
178 else
179 {
180 field1.init(varList1[varID]);
181 hset.getVarLevelPercentiles(field1, varID, levelID, pn);
182 cdo_write_record(streamID4, field1);
183 }
184 }
185
186 if (nrecs == 0) break;
187 otsID++;
188 }
189
190 cdo_stream_close(streamID4);
191 cdo_stream_close(streamID3);
192 cdo_stream_close(streamID2);
193 cdo_stream_close(streamID1);
194 }
195
196 void *
Timpctl(void * process)197 Timpctl(void *process)
198 {
199 cdo_initialize(process);
200
201 // clang-format off
202 cdo_operator_add("timpctl", FieldFunc_Pctl, 31, nullptr);
203 cdo_operator_add("yearpctl", FieldFunc_Pctl, 10, nullptr);
204 cdo_operator_add("monpctl", FieldFunc_Pctl, 8, nullptr);
205 cdo_operator_add("daypctl", FieldFunc_Pctl, 6, nullptr);
206 cdo_operator_add("hourpctl", FieldFunc_Pctl, 4, nullptr);
207 // clang-format on
208
209 timpctl(cdo_operator_id());
210
211 cdo_finish();
212
213 return nullptr;
214 }
215