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 Ydaypctl ydaypctl Multi-year daily percentiles
15 */
16
17 #include <cdi.h>
18
19 #include "cdo_options.h"
20 #include "cdo_vlist.h"
21 #include "datetime.h"
22 #include "process_int.h"
23 #include "param_conversion.h"
24 #include "percentiles_hist.h"
25
26 void *
Ydaypctl(void * process)27 Ydaypctl(void *process)
28 {
29 constexpr int MaxDays = 373;
30 int64_t vdates1[MaxDays] = { 0 }, vdates2[MaxDays] = { 0 };
31 int vtimes1[MaxDays] = { 0 };
32 long nsets[MaxDays] = { 0 };
33 std::vector<bool> vars1(MaxDays, false);
34 HistogramSet hsets[MaxDays];
35
36 cdo_initialize(process);
37 cdo_operator_add("ydaypctl", FieldFunc_Pctl, 0, nullptr);
38
39 operator_input_arg("percentile number");
40 const auto pn = parameter_to_double(cdo_operator_argv(0));
41
42 const auto streamID1 = cdo_open_read(0);
43 const auto streamID2 = cdo_open_read(1);
44 const auto streamID3 = cdo_open_read(2);
45
46 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
47 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
48 const auto vlistID3 = cdo_stream_inq_vlist(streamID3);
49 const auto vlistID4 = vlistDuplicate(vlistID1);
50
51 vlist_compare(vlistID1, vlistID2, CMP_ALL);
52 vlist_compare(vlistID1, vlistID3, CMP_ALL);
53
54 const auto taxisID1 = vlistInqTaxis(vlistID1);
55 const auto taxisID2 = vlistInqTaxis(vlistID2);
56 const auto taxisID3 = vlistInqTaxis(vlistID3);
57 // TODO - check that time axes 2 and 3 are equal
58
59 const auto taxisID4 = taxisDuplicate(taxisID1);
60 if (taxisHasBounds(taxisID4)) taxisDeleteBounds(taxisID4);
61 vlistDefTaxis(vlistID4, taxisID4);
62
63 const auto streamID4 = cdo_open_write(3);
64 cdo_def_vlist(streamID4, vlistID4);
65
66 const auto ntsteps = vlistNtsteps(vlistID1);
67 const auto nvars = vlistNvars(vlistID1);
68
69 const auto maxrecs = vlistNrecs(vlistID1);
70 std::vector<RecordInfo> recList(maxrecs);
71
72 FieldVector constFields(maxrecs);
73
74 Field field1, field2;
75
76 VarList varList1;
77 varListInit(varList1, vlistID1);
78
79 int tsID = 0;
80 while (true)
81 {
82 const auto nrecs = cdo_stream_inq_timestep(streamID2, tsID);
83 if (nrecs == 0) break;
84
85 if (nrecs != cdo_stream_inq_timestep(streamID3, tsID))
86 cdo_abort("Number of records at time step %d of %s and %s differ!", tsID + 1, cdo_get_stream_name(1),
87 cdo_get_stream_name(2));
88
89 const auto vdate = taxisInqVdate(taxisID2);
90 // const auto vtime = taxisInqVtime(taxisID2);
91
92 if (vdate != taxisInqVdate(taxisID3))
93 cdo_abort("Verification dates at time step %d of %s and %s differ!", tsID + 1, cdo_get_stream_name(1),
94 cdo_get_stream_name(2));
95
96 // if (Options::cdoVerbose) cdo_print("process timestep: %d %d %d", tsID + 1, vdate, vtime);
97
98 const auto dayoy = decode_day_of_year(vdate);
99 if (dayoy < 0 || dayoy >= MaxDays) cdo_abort("Day %d out of range!", dayoy);
100
101 vdates2[dayoy] = vdate;
102
103 if (!vars1[dayoy])
104 {
105 vars1[dayoy] = true;
106 hsets[dayoy].create(nvars, ntsteps);
107 for (int varID = 0; varID < nvars; varID++)
108 hsets[dayoy].createVarLevels(varID, varList1[varID].nlevels, varList1[varID].gridsize);
109 }
110
111 for (int recID = 0; recID < nrecs; recID++)
112 {
113 int varID, levelID;
114 cdo_inq_record(streamID2, &varID, &levelID);
115 field1.init(varList1[varID]);
116 cdo_read_record(streamID2, field1);
117
118 cdo_inq_record(streamID3, &varID, &levelID);
119 field2.init(varList1[varID]);
120 cdo_read_record(streamID3, field2);
121
122 hsets[dayoy].defVarLevelBounds(varID, levelID, field1, field2);
123 }
124
125 tsID++;
126 }
127
128 tsID = 0;
129 while (true)
130 {
131 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
132 if (nrecs == 0) break;
133
134 const auto vdate = taxisInqVdate(taxisID1);
135 const auto vtime = taxisInqVtime(taxisID1);
136
137 // if (Options::cdoVerbose) cdo_print("process timestep: %d %d %d", tsID + 1, vdate, vtime);
138
139 const auto dayoy = decode_day_of_year(vdate);
140 if (dayoy < 0 || dayoy >= MaxDays) cdo_abort("Day %d out of range!", dayoy);
141
142 vdates1[dayoy] = vdate;
143 vtimes1[dayoy] = vtime;
144
145 if (!vars1[dayoy]) cdo_abort("No data for day %d in %s and %s", dayoy, cdo_get_stream_name(1), cdo_get_stream_name(2));
146
147 for (int recID = 0; recID < nrecs; recID++)
148 {
149 int varID, levelID;
150 cdo_inq_record(streamID1, &varID, &levelID);
151
152 if (tsID == 0)
153 {
154 recList[recID].varID = varID;
155 recList[recID].levelID = levelID;
156 recList[recID].lconst = varList1[varID].timetype == TIME_CONSTANT;
157 }
158
159 if (tsID == 0 && recList[recID].lconst)
160 {
161 constFields[recID].init(varList1[varID]);
162 cdo_read_record(streamID1, constFields[recID]);
163 }
164 else
165 {
166 field1.init(varList1[varID]);
167 cdo_read_record(streamID1, field1);
168
169 hsets[dayoy].addVarLevelValues(varID, levelID, field1);
170 }
171 }
172
173 nsets[dayoy]++;
174 tsID++;
175 }
176
177 int otsID = 0;
178 for (int dayoy = 0; dayoy < MaxDays; dayoy++)
179 if (nsets[dayoy])
180 {
181 if (decode_month_and_day(vdates1[dayoy]) != decode_month_and_day(vdates2[dayoy]))
182 cdo_abort("Verification dates for the day %d of %s and %s are different!", dayoy, cdo_get_stream_name(0),
183 cdo_get_stream_name(1));
184
185 taxisDefVdate(taxisID4, vdates1[dayoy]);
186 taxisDefVtime(taxisID4, vtimes1[dayoy]);
187 cdo_def_timestep(streamID4, otsID);
188
189 for (int recID = 0; recID < maxrecs; recID++)
190 {
191 if (otsID && recList[recID].lconst) continue;
192
193 const auto varID = recList[recID].varID;
194 const auto levelID = recList[recID].levelID;
195 cdo_def_record(streamID4, varID, levelID);
196
197 if (recList[recID].lconst)
198 {
199 cdo_write_record(streamID4, constFields[recID]);
200 }
201 else
202 {
203 field1.init(varList1[varID]);
204 hsets[dayoy].getVarLevelPercentiles(field1, varID, levelID, pn);
205 cdo_write_record(streamID4, field1);
206 }
207 }
208
209 otsID++;
210 }
211
212 cdo_stream_close(streamID4);
213 cdo_stream_close(streamID3);
214 cdo_stream_close(streamID2);
215 cdo_stream_close(streamID1);
216
217 cdo_finish();
218
219 return nullptr;
220 }
221