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
6 */
7
8 /*
9 This module contains the following operators:
10
11 Ymonstat ymonrange Multi-year monthly range
12 Ymonstat ymonmin Multi-year monthly minimum
13 Ymonstat ymonmax Multi-year monthly maximum
14 Ymonstat ymonsum Multi-year monthly sum
15 Ymonstat ymonmean Multi-year monthly mean
16 Ymonstat ymonavg Multi-year monthly average
17 Ymonstat ymonvar Multi-year monthly variance
18 Ymonstat ymonvar1 Multi-year monthly variance [Normalize by (n-1)]
19 Ymonstat ymonstd Multi-year monthly standard deviation
20 Ymonstat ymonstd1 Multi-year monthly standard deviation [Normalize by (n-1)]
21 */
22
23 #include <cdi.h>
24
25 #include "cdo_options.h"
26 #include "datetime.h"
27 #include "process_int.h"
28
29 static void
addOperators(void)30 addOperators(void)
31 {
32 // clang-format off
33 cdo_operator_add("ymonrange", FieldFunc_Range, 0, nullptr);
34 cdo_operator_add("ymonmin", FieldFunc_Min, 0, nullptr);
35 cdo_operator_add("ymonmax", FieldFunc_Max, 0, nullptr);
36 cdo_operator_add("ymonsum", FieldFunc_Sum, 0, nullptr);
37 cdo_operator_add("ymonmean", FieldFunc_Mean, 0, nullptr);
38 cdo_operator_add("ymonavg", FieldFunc_Avg, 0, nullptr);
39 cdo_operator_add("ymonvar", FieldFunc_Var, 0, nullptr);
40 cdo_operator_add("ymonvar1", FieldFunc_Var1, 0, nullptr);
41 cdo_operator_add("ymonstd", FieldFunc_Std, 0, nullptr);
42 cdo_operator_add("ymonstd1", FieldFunc_Std1, 0, nullptr);
43 // clang-format on
44 }
45
46 void *
Ymonstat(void * process)47 Ymonstat(void *process)
48 {
49 const TimeStat timestat_date = TimeStat::LAST;
50 constexpr int MaxMonths = 17;
51 int month_nsets[MaxMonths] = { 0 };
52 int mon[MaxMonths] = { 0 };
53 int nmon = 0;
54 FieldVector2D vars1[MaxMonths], vars2[MaxMonths], samp1[MaxMonths];
55
56 cdo_initialize(process);
57
58 addOperators();
59
60 const auto operatorID = cdo_operator_id();
61 const auto operfunc = cdo_operator_f1(operatorID);
62
63 const auto lrange = (operfunc == FieldFunc_Range);
64 const auto lmean = (operfunc == FieldFunc_Mean || operfunc == FieldFunc_Avg);
65 const auto lstd = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Std1);
66 const auto lvarstd = (lstd || operfunc == FieldFunc_Var || operfunc == FieldFunc_Var1);
67 const auto lvars2 = (lvarstd || lrange);
68 const int divisor = (operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
69
70 auto field2_stdvar_func = lstd ? field2_std : field2_var;
71 auto fieldc_stdvar_func = lstd ? fieldc_std : fieldc_var;
72
73 operator_check_argc(0);
74
75 const auto streamID1 = cdo_open_read(0);
76
77 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
78 const auto vlistID2 = vlistDuplicate(vlistID1);
79
80 const auto taxisID1 = vlistInqTaxis(vlistID1);
81 const auto taxisID2 = taxisDuplicate(taxisID1);
82 taxisWithBounds(taxisID2);
83 if (taxisInqType(taxisID2) == TAXIS_FORECAST) taxisDefType(taxisID2, TAXIS_RELATIVE);
84 vlistDefTaxis(vlistID2, taxisID2);
85
86 const auto streamID2 = cdo_open_write(1);
87 cdo_def_vlist(streamID2, vlistID2);
88
89 const auto maxrecs = vlistNrecs(vlistID1);
90 std::vector<RecordInfo> recList(maxrecs);
91
92 DateTimeList dtlists[MaxMonths];
93 for (int month = 0; month < MaxMonths; ++month)
94 {
95 dtlists[month].set_stat(timestat_date);
96 dtlists[month].set_calendar(taxisInqCalendar(taxisID1));
97 }
98
99 VarList varList;
100 varListInit(varList, vlistID1);
101
102 int VARS_MEMTYPE = 0;
103 if ((operfunc == FieldFunc_Min) || (operfunc == FieldFunc_Max)) VARS_MEMTYPE = FIELD_NAT;
104
105 Field field;
106
107 int tsID = 0;
108 int otsID = 0;
109 while (true)
110 {
111 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
112 if (nrecs == 0) break;
113
114 const auto vdate = taxisInqVdate(taxisID1);
115 const auto vtime = taxisInqVtime(taxisID1);
116
117 if (Options::cdoVerbose) cdo_print("process timestep: %d %d %d", tsID + 1, vdate, vtime);
118
119 const auto month = decode_month(vdate);
120 if (month < 0 || month >= MaxMonths) cdo_abort("Month %d out of range!", month);
121
122 dtlists[month].taxis_set_next_timestep(taxisID1);
123
124 if (!vars1[month].size())
125 {
126 mon[nmon++] = month;
127 fields_from_vlist(vlistID1, samp1[month]);
128 fields_from_vlist(vlistID1, vars1[month], FIELD_VEC | VARS_MEMTYPE);
129 if (lvars2) fields_from_vlist(vlistID1, vars2[month], FIELD_VEC);
130 }
131
132 for (int recID = 0; recID < nrecs; recID++)
133 {
134 int varID, levelID;
135 cdo_inq_record(streamID1, &varID, &levelID);
136
137 if (tsID == 0)
138 {
139 recList[recID].varID = varID;
140 recList[recID].levelID = levelID;
141 recList[recID].lconst = (varList[varID].timetype == TIME_CONSTANT);
142 }
143
144 auto &rsamp1 = samp1[month][varID][levelID];
145 auto &rvars1 = vars1[month][varID][levelID];
146
147 const auto nsets = month_nsets[month];
148
149 if (nsets == 0)
150 {
151 cdo_read_record(streamID1, rvars1);
152 if (lrange)
153 {
154 vars2[month][varID][levelID].nmiss = rvars1.nmiss;
155 vars2[month][varID][levelID].vec_d = rvars1.vec_d;
156 }
157
158 if (rvars1.nmiss || !rsamp1.empty())
159 {
160 if (rsamp1.empty()) rsamp1.resize(rvars1.size);
161 field2_vinit(rsamp1, rvars1);
162 }
163 }
164 else
165 {
166 field.init(varList[varID]);
167 cdo_read_record(streamID1, field);
168
169 if (field.nmiss || !rsamp1.empty())
170 {
171 if (rsamp1.empty()) rsamp1.resize(rvars1.size, nsets);
172 field2_vincr(rsamp1, field);
173 }
174
175 // clang-format off
176 if (lvarstd) field2_sumsumq(rvars1, vars2[month][varID][levelID], field);
177 else if (lrange) field2_maxmin(rvars1, vars2[month][varID][levelID], field);
178 else field2_function(rvars1, field, operfunc);
179 // clang-format on
180 }
181 }
182
183 if (month_nsets[month] == 0 && lvarstd)
184 for (int recID = 0; recID < maxrecs; recID++)
185 {
186 if (recList[recID].lconst) continue;
187
188 const auto varID = recList[recID].varID;
189 const auto levelID = recList[recID].levelID;
190 field2_moq(vars2[month][varID][levelID], vars1[month][varID][levelID]);
191 }
192
193 month_nsets[month]++;
194 tsID++;
195 }
196
197 if (nmon == 12)
198 {
199 int smon = 0;
200 for (int month = 1; month <= 12; month++)
201 if (month_nsets[month]) smon++;
202 if (smon == 12)
203 for (int month = 1; month <= 12; month++) mon[month - 1] = month;
204 }
205
206 for (int i = 0; i < nmon; i++)
207 {
208 const auto month = mon[i];
209 const auto nsets = month_nsets[month];
210 if (nsets == 0) cdo_abort("Internal problem, nsets[%d] not defined!", month);
211
212 for (int recID = 0; recID < maxrecs; recID++)
213 {
214 if (recList[recID].lconst) continue;
215
216 const auto varID = recList[recID].varID;
217 const auto levelID = recList[recID].levelID;
218 const auto &rsamp1 = samp1[month][varID][levelID];
219 auto &rvars1 = vars1[month][varID][levelID];
220
221 if (lmean)
222 {
223 if (!rsamp1.empty())
224 field2_div(rvars1, rsamp1);
225 else
226 fieldc_div(rvars1, (double) nsets);
227 }
228 else if (lvarstd)
229 {
230 if (!rsamp1.empty())
231 field2_stdvar_func(rvars1, vars2[month][varID][levelID], rsamp1, divisor);
232 else
233 fieldc_stdvar_func(rvars1, vars2[month][varID][levelID], nsets, divisor);
234 }
235 else if (lrange)
236 {
237 field2_sub(rvars1, vars2[month][varID][levelID]);
238 }
239 }
240
241 dtlists[month].stat_taxis_def_timestep(taxisID2);
242 cdo_def_timestep(streamID2, otsID);
243
244 for (int recID = 0; recID < maxrecs; recID++)
245 {
246 if (otsID && recList[recID].lconst) continue;
247
248 const auto varID = recList[recID].varID;
249 const auto levelID = recList[recID].levelID;
250 auto &rvars1 = vars1[month][varID][levelID];
251
252 cdo_def_record(streamID2, varID, levelID);
253 cdo_write_record(streamID2, rvars1);
254 }
255
256 otsID++;
257 }
258
259 cdo_stream_close(streamID2);
260 cdo_stream_close(streamID1);
261
262 cdo_finish();
263
264 return nullptr;
265 }
266