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 #include <cdi.h>
9
10 #include "process_int.h"
11
12 static void
checkUniqueZaxis(int vlistID)13 checkUniqueZaxis(int vlistID)
14 {
15 const auto nzaxis = vlistNzaxis(vlistID);
16 const auto zaxisID = vlistZaxis(vlistID, 0);
17 const auto nlevels = zaxisInqSize(zaxisID);
18 for (int index = 1; index < nzaxis; ++index)
19 {
20 if (nlevels != zaxisInqSize(vlistZaxis(vlistID, index))) cdo_abort("Number of level differ!");
21 }
22 }
23
24 static void
checkUniqueGridsize(int vlistID)25 checkUniqueGridsize(int vlistID)
26 {
27 const auto ngrids = vlistNgrids(vlistID);
28 const auto gridID = vlistGrid(vlistID, 0);
29 const auto gridsize = gridInqSize(gridID);
30 for (int index = 0; index < ngrids; index++)
31 {
32 if (gridsize != gridInqSize(vlistGrid(vlistID, index))) cdo_abort("Horizontal gridsize differ!");
33 }
34 }
35
36 static void
setAttributes(const VarList & varList1,int vlistID1,int vlistID2,int varID2,int operatorID)37 setAttributes(const VarList &varList1, int vlistID1, int vlistID2, int varID2, int operatorID)
38 {
39 auto paramIsEqual = true;
40 char name[CDI_MAX_NAME], zname[CDI_MAX_NAME];
41 strcpy(name, varList1[0].name);
42 const auto param = varList1[0].param;
43 const auto nvars = vlistNvars(vlistID1);
44 for (int varID = 1; varID < nvars; ++varID)
45 {
46 strcpy(zname, varList1[varID].name);
47 const auto zparam = varList1[varID].param;
48 if (param != zparam || !cdo_cmpstr(name, zname))
49 {
50 paramIsEqual = false;
51 break;
52 }
53 }
54
55 if (!paramIsEqual) strcpy(name, cdo_operator_name(operatorID));
56 cdiDefKeyString(vlistID2, varID2, CDI_KEY_NAME, name);
57 if (paramIsEqual)
58 {
59 if (param >= 0) vlistDefVarParam(vlistID2, varID2, param);
60 zname[0] = 0;
61 vlistInqVarLongname(vlistID1, 0, zname);
62 if (*zname) cdiDefKeyString(vlistID2, varID2, CDI_KEY_LONGNAME, zname);
63 zname[0] = 0;
64 vlistInqVarUnits(vlistID1, 0, zname);
65 if (*zname) cdiDefKeyString(vlistID2, varID2, CDI_KEY_UNITS, zname);
66 }
67 }
68
69 static void
addOperators(void)70 addOperators(void)
71 {
72 // clang-format off
73 cdo_operator_add("varsrange", FieldFunc_Range, 0, nullptr);
74 cdo_operator_add("varsmin", FieldFunc_Min, 0, nullptr);
75 cdo_operator_add("varsmax", FieldFunc_Max, 0, nullptr);
76 cdo_operator_add("varssum", FieldFunc_Sum, 0, nullptr);
77 cdo_operator_add("varsmean", FieldFunc_Mean, 0, nullptr);
78 cdo_operator_add("varsavg", FieldFunc_Avg, 0, nullptr);
79 cdo_operator_add("varsvar", FieldFunc_Var, 0, nullptr);
80 cdo_operator_add("varsvar1", FieldFunc_Var1, 0, nullptr);
81 cdo_operator_add("varsstd", FieldFunc_Std, 0, nullptr);
82 cdo_operator_add("varsstd1", FieldFunc_Std1, 0, nullptr);
83 // clang-format on
84 }
85
86 void *
Varsstat(void * process)87 Varsstat(void *process)
88 {
89 cdo_initialize(process);
90
91 addOperators();
92
93 const auto operatorID = cdo_operator_id();
94 const auto operfunc = cdo_operator_f1(operatorID);
95
96 operator_check_argc(0);
97
98 const auto lrange = (operfunc == FieldFunc_Range);
99 const auto lmean = (operfunc == FieldFunc_Mean || operfunc == FieldFunc_Avg);
100 const auto lstd = (operfunc == FieldFunc_Std || operfunc == FieldFunc_Std1);
101 const auto lvarstd = (lstd || operfunc == FieldFunc_Var || operfunc == FieldFunc_Var1);
102 const auto lvars2 = (lvarstd || lrange);
103 const int divisor = (operfunc == FieldFunc_Std1 || operfunc == FieldFunc_Var1);
104
105 auto field2_stdvar_func = lstd ? field2_std : field2_var;
106 auto fieldc_stdvar_func = lstd ? fieldc_std : fieldc_var;
107
108 const auto streamID1 = cdo_open_read(0);
109 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
110
111 VarList varList1;
112 varListInit(varList1, vlistID1);
113
114 checkUniqueZaxis(vlistID1);
115 const auto zaxisID = vlistZaxis(vlistID1, 0);
116 const auto nlevels = zaxisInqSize(zaxisID);
117
118 checkUniqueGridsize(vlistID1);
119 const auto gridID = vlistGrid(vlistID1, 0);
120 const auto gridsize = gridInqSize(gridID);
121
122 const auto timetype = varList1[0].timetype;
123 const auto nvars = vlistNvars(vlistID1);
124 for (int varID = 1; varID < nvars; ++varID)
125 {
126 if (timetype != varList1[varID].timetype) cdo_abort("Number of timesteps differ!");
127 }
128
129 const auto vlistID2 = vlistCreate();
130 vlistDefNtsteps(vlistID2, vlistNtsteps(vlistID1));
131
132 const auto varID2 = vlistDefVar(vlistID2, gridID, zaxisID, timetype);
133 setAttributes(varList1, vlistID1, vlistID2, varID2, operatorID);
134
135 const auto taxisID1 = vlistInqTaxis(vlistID1);
136 const auto taxisID2 = taxisDuplicate(taxisID1);
137 vlistDefTaxis(vlistID2, taxisID2);
138
139 const auto streamID2 = cdo_open_write(1);
140 cdo_def_vlist(streamID2, vlistID2);
141
142 Field field;
143
144 FieldVector vars1(nlevels), samp1(nlevels), vars2;
145 if (lvars2) vars2.resize(nlevels);
146
147 for (int levelID = 0; levelID < nlevels; levelID++)
148 {
149 const auto missval = varList1[0].missval;
150
151 samp1[levelID].grid = gridID;
152 samp1[levelID].missval = missval;
153 samp1[levelID].memType = MemType::Double;
154 vars1[levelID].grid = gridID;
155 vars1[levelID].missval = missval;
156 vars1[levelID].memType = MemType::Double;
157 vars1[levelID].resize(gridsize);
158 if (lvars2)
159 {
160 vars2[levelID].grid = gridID;
161 vars2[levelID].missval = missval;
162 vars2[levelID].memType = MemType::Double;
163 vars2[levelID].resize(gridsize);
164 }
165 }
166
167 int tsID = 0;
168 while (true)
169 {
170 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
171 if (nrecs == 0) break;
172
173 cdo_taxis_copy_timestep(taxisID2, taxisID1);
174 cdo_def_timestep(streamID2, tsID);
175
176 for (int recID = 0; recID < nrecs; recID++)
177 {
178 int varID, levelID;
179 cdo_inq_record(streamID1, &varID, &levelID);
180
181 auto &rsamp1 = samp1[levelID];
182 auto &rvars1 = vars1[levelID];
183
184 rvars1.nsamp++;
185 if (lrange) vars2[levelID].nsamp++;
186
187 if (varID == 0)
188 {
189 cdo_read_record(streamID1, rvars1);
190 if (lrange)
191 {
192 vars2[levelID].nmiss = rvars1.nmiss;
193 vars2[levelID].vec_d = rvars1.vec_d;
194 }
195
196 if (lvarstd) field2_moq(vars2[levelID], rvars1);
197
198 if (rvars1.nmiss || !rsamp1.empty())
199 {
200 if (rsamp1.empty()) rsamp1.resize(rvars1.size);
201 field2_vinit(rsamp1, rvars1);
202 }
203 }
204 else
205 {
206 field.init(varList1[varID]);
207 cdo_read_record(streamID1, field);
208
209 if (field.nmiss || !rsamp1.empty())
210 {
211 if (rsamp1.empty()) rsamp1.resize(rvars1.size, rvars1.nsamp);
212 field2_vincr(rsamp1, field);
213 }
214
215 // clang-format off
216 if (lvarstd) field2_sumsumq(rvars1, vars2[levelID], field);
217 else if (lrange) field2_maxmin(rvars1, vars2[levelID], field);
218 else field2_function(rvars1, field, operfunc);
219 // clang-format on
220 }
221 }
222
223 for (int levelID = 0; levelID < nlevels; levelID++)
224 {
225 const auto &rsamp1 = samp1[levelID];
226 auto &rvars1 = vars1[levelID];
227
228 if (rvars1.nsamp)
229 {
230 if (lmean)
231 {
232 if (!rsamp1.empty())
233 field2_div(rvars1, rsamp1);
234 else
235 fieldc_div(rvars1, (double) rvars1.nsamp);
236 }
237 else if (lvarstd)
238 {
239 if (!rsamp1.empty())
240 field2_stdvar_func(rvars1, vars2[levelID], rsamp1, divisor);
241 else
242 fieldc_stdvar_func(rvars1, vars2[levelID], rsamp1.nsamp, divisor);
243 }
244 else if (lrange)
245 {
246 field2_sub(rvars1, vars2[levelID]);
247 }
248
249 cdo_def_record(streamID2, 0, levelID);
250 cdo_write_record(streamID2, rvars1);
251 rvars1.nsamp = 0;
252 }
253 }
254
255 tsID++;
256 }
257
258 cdo_stream_close(streamID2);
259 cdo_stream_close(streamID1);
260
261 vlistDestroy(vlistID2);
262
263 cdo_finish();
264
265 return nullptr;
266 }
267