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 Ralf Quast
6
7 */
8
9 /*
10 This module contains the following operators:
11
12 Runpctl runpctl Running percentiles
13 */
14
15 #include <cdi.h>
16
17 #include "process_int.h"
18 #include "param_conversion.h"
19 #include "percentiles.h"
20 #include "datetime.h"
21
22 void *
Runpctl(void * process)23 Runpctl(void *process)
24 {
25 const auto timestat_date = TimeStat::MEAN;
26
27 cdo_initialize(process);
28
29 operator_input_arg("percentile number, number of timesteps");
30 operator_check_argc(2);
31 const auto pn = parameter_to_double(cdo_operator_argv(0));
32 const auto ndates = parameter_to_int(cdo_operator_argv(1));
33
34 const auto streamID1 = cdo_open_read(0);
35
36 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
37 const auto vlistID2 = vlistDuplicate(vlistID1);
38
39 const auto taxisID1 = vlistInqTaxis(vlistID1);
40 const auto taxisID2 = taxisDuplicate(taxisID1);
41 taxisWithBounds(taxisID2);
42 vlistDefTaxis(vlistID2, taxisID2);
43
44 const auto streamID2 = cdo_open_write(1);
45 cdo_def_vlist(streamID2, vlistID2);
46
47 const auto nvars = vlistNvars(vlistID1);
48
49 const auto maxrecs = vlistNrecs(vlistID1);
50 std::vector<RecordInfo> recList(maxrecs);
51
52 DateTimeList dtlist;
53 dtlist.set_stat(timestat_date);
54 dtlist.set_calendar(taxisInqCalendar(taxisID1));
55
56 Varray<float> array_f(ndates);
57 Varray<double> array_d(ndates);
58 FieldVector3D vars1(ndates + 1);
59
60 for (int its = 0; its < ndates; its++) fields_from_vlist(vlistID1, vars1[its]);
61
62 VarList varList1;
63 varListInit(varList1, vlistID1);
64
65 int tsID;
66 for (tsID = 0; tsID < ndates; tsID++)
67 {
68 auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
69 if (nrecs == 0) cdo_abort("File has less than %d timesteps!", ndates);
70
71 dtlist.taxis_inq_timestep(taxisID1, tsID);
72
73 for (int recID = 0; recID < nrecs; recID++)
74 {
75 int varID, levelID;
76 cdo_inq_record(streamID1, &varID, &levelID);
77
78 if (tsID == 0)
79 {
80 recList[recID].varID = varID;
81 recList[recID].levelID = levelID;
82 recList[recID].lconst = (varList1[varID].timetype == TIME_CONSTANT);
83 }
84
85 auto &field = vars1[tsID][varID][levelID];
86 field.init(varList1[varID]);
87 cdo_read_record(streamID1, field);
88 }
89 }
90
91 int otsID = 0;
92 while (true)
93 {
94 for (int varID = 0; varID < nvars; varID++)
95 {
96 if (varList1[varID].timetype == TIME_CONSTANT) continue;
97
98 const auto gridsize = varList1[varID].gridsize;
99 const auto nlevels = varList1[varID].nlevels;
100 const auto missval = varList1[varID].missval;
101
102 for (int levelID = 0; levelID < nlevels; levelID++)
103 {
104 auto &field1 = vars1[0][varID][levelID];
105 size_t nmiss = 0;
106 if (field1.memType == MemType::Float)
107 {
108 for (size_t i = 0; i < gridsize; i++)
109 {
110 int j = 0;
111 for (int inp = 0; inp < ndates; inp++)
112 {
113 const auto val = vars1[inp][varID][levelID].vec_f[i];
114 if (!DBL_IS_EQUAL(val, (float) missval)) array_f[j++] = val;
115 }
116
117 if (j > 0)
118 {
119 field1.vec_f[i] = percentile(array_f.data(), j, pn);
120 }
121 else
122 {
123 field1.vec_f[i] = missval;
124 nmiss++;
125 }
126 }
127 }
128 else
129 {
130 for (size_t i = 0; i < gridsize; i++)
131 {
132 int j = 0;
133 for (int inp = 0; inp < ndates; inp++)
134 {
135 const auto val = vars1[inp][varID][levelID].vec_d[i];
136 if (!DBL_IS_EQUAL(val, missval)) array_d[j++] = val;
137 }
138
139 if (j > 0)
140 {
141 field1.vec_d[i] = percentile(array_d.data(), j, pn);
142 }
143 else
144 {
145 field1.vec_d[i] = missval;
146 nmiss++;
147 }
148 }
149 }
150 field1.nmiss = nmiss;
151 }
152 }
153
154 dtlist.stat_taxis_def_timestep(taxisID2, ndates);
155 cdo_def_timestep(streamID2, otsID);
156
157 for (int recID = 0; recID < maxrecs; recID++)
158 {
159 if (otsID && recList[recID].lconst) continue;
160
161 const auto varID = recList[recID].varID;
162 const auto levelID = recList[recID].levelID;
163
164 cdo_def_record(streamID2, varID, levelID);
165 auto &field1 = vars1[0][varID][levelID];
166 cdo_write_record(streamID2, field1);
167 }
168
169 otsID++;
170
171 dtlist.shift();
172
173 vars1[ndates] = vars1[0];
174 for (int inp = 0; inp < ndates; inp++) vars1[inp] = vars1[inp + 1];
175
176 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
177 if (nrecs == 0) break;
178
179 dtlist.taxis_inq_timestep(taxisID1, ndates - 1);
180
181 for (int recID = 0; recID < nrecs; recID++)
182 {
183 int varID, levelID;
184 cdo_inq_record(streamID1, &varID, &levelID);
185 auto &fieldN = vars1[ndates - 1][varID][levelID];
186 cdo_read_record(streamID1, fieldN);
187 }
188
189 tsID++;
190 }
191
192 cdo_stream_close(streamID2);
193 cdo_stream_close(streamID1);
194
195 cdo_finish();
196
197 return nullptr;
198 }
199