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 Invert invertlat Invert latitude
12 Invert invertlon Invert longitude
13 Invert invertlatdes Invert latitude description
14 Invert invertlondes Invert longitude description
15 Invert invertlatdata Invert latitude data
16 Invert invertlondata Invert longitude data
17 */
18
19 #include <cdi.h>
20
21 #include "process_int.h"
22
23 static void
invertLonDes(int vlistID)24 invertLonDes(int vlistID)
25 {
26 const auto ngrids = vlistNgrids(vlistID);
27 for (int index = 0; index < ngrids; index++)
28 {
29 const auto gridID1 = vlistGrid(vlistID, index);
30 const auto gridID2 = gridDuplicate(gridID1);
31
32 const auto gridtype = gridInqType(gridID1);
33
34 if (!(gridtype == GRID_GENERIC || gridtype == GRID_GAUSSIAN || gridtype == GRID_PROJECTION || gridtype == GRID_LONLAT
35 || gridtype == GRID_CURVILINEAR))
36 cdo_abort("Unsupported gridtype: %s!", gridNamePtr(gridtype));
37
38 if (gridInqXvals(gridID1, nullptr))
39 {
40 const auto nlon = gridInqXsize(gridID1);
41 const auto nlat = gridInqYsize(gridID1);
42 const auto size = (gridtype == GRID_CURVILINEAR) ? nlon * nlat : nlon;
43
44 Varray<double> xv1(size), xv2(size);
45
46 gridInqXvals(gridID1, xv1.data());
47
48 if (gridtype == GRID_CURVILINEAR)
49 {
50 for (size_t ilat = 0; ilat < nlat; ilat++)
51 for (size_t ilon = 0; ilon < nlon; ilon++) xv2[ilat * nlon + nlon - ilon - 1] = xv1[ilat * nlon + ilon];
52 }
53 else
54 {
55 for (size_t ilon = 0; ilon < nlon; ilon++) xv2[nlon - ilon - 1] = xv1[ilon];
56 }
57
58 gridDefXvals(gridID2, xv2.data());
59 }
60
61 if (gridInqXbounds(gridID1, nullptr))
62 {
63 const auto nlon = gridInqXsize(gridID1);
64 const auto nlat = gridInqYsize(gridID1);
65 const auto nv = gridInqNvertex(gridID1);
66 const auto size = (gridtype == GRID_CURVILINEAR) ? nv * nlon * nlat : nv * nlon;
67
68 Varray<double> xb1(size), xb2(size);
69
70 gridInqXbounds(gridID1, xb1.data());
71
72 if (gridtype == GRID_CURVILINEAR)
73 {
74 for (size_t ilat = 0; ilat < nlat; ilat++)
75 for (size_t ilon = 0; ilon < nlon; ilon++)
76 for (int iv = 0; iv < nv; iv++)
77 xb2[ilat * nlon * nv + (nlon - ilon - 1) * nv + iv] = xb1[ilat * nlon * nv + ilon * nv + iv];
78 }
79 else
80 {
81 for (size_t ilon = 0; ilon < nlon; ilon++)
82 {
83 xb2[nlon * 2 - ilon * 2 - 1] = xb1[ilon * 2];
84 xb2[nlon * 2 - ilon * 2 - 2] = xb1[ilon * 2 + 1];
85 }
86 }
87
88 gridDefXbounds(gridID2, xb2.data());
89 }
90
91 vlistChangeGrid(vlistID, gridID1, gridID2);
92 }
93 }
94
95 static void
invertLatCoord(int gridID)96 invertLatCoord(int gridID)
97 {
98 const auto gridtype = gridInqType(gridID);
99
100 if (gridInqYvals(gridID, nullptr))
101 {
102 const auto nlon = gridInqXsize(gridID);
103 const auto nlat = gridInqYsize(gridID);
104 const auto size = (gridtype == GRID_CURVILINEAR) ? nlon * nlat : nlat;
105
106 Varray<double> yv1(size), yv2(size);
107
108 if (gridtype == GRID_CURVILINEAR)
109 {
110 gridInqXvals(gridID, yv1.data());
111
112 for (size_t ilat = 0; ilat < nlat; ilat++)
113 for (size_t ilon = 0; ilon < nlon; ilon++) yv2[(nlat - ilat - 1) * nlon + ilon] = yv1[ilat * nlon + ilon];
114
115 gridDefXvals(gridID, yv2.data());
116
117 gridInqYvals(gridID, yv1.data());
118
119 for (size_t ilat = 0; ilat < nlat; ilat++)
120 for (size_t ilon = 0; ilon < nlon; ilon++) yv2[(nlat - ilat - 1) * nlon + ilon] = yv1[ilat * nlon + ilon];
121
122 gridDefYvals(gridID, yv2.data());
123 }
124 else
125 {
126 gridInqYvals(gridID, yv1.data());
127
128 for (size_t ilat = 0; ilat < nlat; ilat++) yv2[nlat - ilat - 1] = yv1[ilat];
129
130 gridDefYvals(gridID, yv2.data());
131 }
132 }
133
134 if (gridInqYbounds(gridID, nullptr))
135 {
136 const auto nlon = gridInqXsize(gridID);
137 const auto nlat = gridInqYsize(gridID);
138 const auto nv = gridInqNvertex(gridID);
139 const auto size = (gridtype == GRID_CURVILINEAR) ? nv * nlon * nlat : nv * nlat;
140
141 Varray<double> yb1(size), yb2(size);
142
143 gridInqYbounds(gridID, yb1.data());
144
145 if (gridtype == GRID_CURVILINEAR)
146 {
147 for (size_t ilat = 0; ilat < nlat; ilat++)
148 for (size_t ilon = 0; ilon < nlon; ilon++)
149 for (int iv = 0; iv < nv; iv++)
150 yb2[(nlat - ilat - 1) * nlon * nv + ilon * nv + iv] = yb1[ilat * nlon * nv + ilon * nv + iv];
151 }
152 else
153 {
154 for (size_t ilat = 0; ilat < nlat; ilat++)
155 {
156 yb2[nlat * 2 - ilat * 2 - 1] = yb1[ilat * 2];
157 yb2[nlat * 2 - ilat * 2 - 2] = yb1[ilat * 2 + 1];
158 }
159 }
160
161 gridDefYbounds(gridID, yb2.data());
162 }
163 }
164
165 static void
invertLatDes(int vlistID)166 invertLatDes(int vlistID)
167 {
168 const auto ngrids = vlistNgrids(vlistID);
169 for (int index = 0; index < ngrids; index++)
170 {
171 const auto gridID1 = vlistGrid(vlistID, index);
172 const auto gridID2 = gridDuplicate(gridID1);
173
174 const auto gridtype = gridInqType(gridID1);
175
176 if (!(gridtype == GRID_GENERIC || gridtype == GRID_GAUSSIAN || gridtype == GRID_PROJECTION || gridtype == GRID_LONLAT
177 || gridtype == GRID_CURVILINEAR))
178 cdo_abort("Unsupported gridtype: %s!", gridNamePtr(gridtype));
179
180 invertLatCoord(gridID2);
181
182 const auto projID = gridInqProj(gridID2);
183 if (projID != CDI_UNDEFID) invertLatCoord(projID);
184
185 vlistChangeGrid(vlistID, gridID1, gridID2);
186 }
187 }
188
189 static void
invertLonData(double * array1,double * array2,int gridID1)190 invertLonData(double *array1, double *array2, int gridID1)
191 {
192 const auto nlon = gridInqXsize(gridID1);
193 const auto nlat = gridInqYsize(gridID1);
194
195 if (nlat > 0)
196 {
197 std::vector<double *> field1(nlat), field2(nlat);
198
199 for (size_t ilat = 0; ilat < nlat; ilat++)
200 {
201 field1[ilat] = array1 + ilat * nlon;
202 field2[ilat] = array2 + ilat * nlon;
203 }
204
205 for (size_t ilat = 0; ilat < nlat; ilat++)
206 for (size_t ilon = 0; ilon < nlon; ilon++) field2[ilat][nlon - ilon - 1] = field1[ilat][ilon];
207 }
208 else
209 {
210 array2[0] = array1[0];
211 }
212 }
213
214 static void
invertLatData(double * array1,double * array2,int gridID1)215 invertLatData(double *array1, double *array2, int gridID1)
216 {
217 const auto nlon = gridInqXsize(gridID1);
218 const auto nlat = gridInqYsize(gridID1);
219
220 if (nlat > 0)
221 {
222 std::vector<double *> field1(nlat), field2(nlat);
223
224 for (size_t ilat = 0; ilat < nlat; ilat++)
225 {
226 field1[ilat] = array1 + ilat * nlon;
227 field2[ilat] = array2 + ilat * nlon;
228 }
229
230 for (size_t ilat = 0; ilat < nlat; ilat++) array_copy(nlon, field1[ilat], field2[nlat - ilat - 1]);
231 }
232 else
233 {
234 array2[0] = array1[0];
235 }
236 }
237
238 void *
Invert(void * process)239 Invert(void *process)
240 {
241 enum
242 {
243 func_fld,
244 func_all,
245 func_hrd,
246 func_lon,
247 func_lat
248 };
249
250 cdo_initialize(process);
251
252 // clang-format off
253 cdo_operator_add("invertlat", func_all, func_lat, nullptr);
254 cdo_operator_add("invertlon", func_all, func_lon, nullptr);
255 cdo_operator_add("invertlatdes", func_hrd, func_lat, nullptr);
256 cdo_operator_add("invertlondes", func_hrd, func_lon, nullptr);
257 cdo_operator_add("invertlatdata", func_fld, func_lat, nullptr);
258 cdo_operator_add("invertlondata", func_fld, func_lon, nullptr);
259 // clang-format on
260
261 const auto operatorID = cdo_operator_id();
262 const auto operfunc1 = cdo_operator_f1(operatorID);
263 const auto operfunc2 = cdo_operator_f2(operatorID);
264
265 operator_check_argc(0);
266
267 const auto streamID1 = cdo_open_read(0);
268
269 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
270 const auto vlistID2 = vlistDuplicate(vlistID1);
271
272 const auto taxisID1 = vlistInqTaxis(vlistID1);
273 const auto taxisID2 = taxisDuplicate(taxisID1);
274 vlistDefTaxis(vlistID2, taxisID2);
275
276 if (operfunc1 == func_all || operfunc1 == func_hrd)
277 {
278 if (operfunc2 == func_lat)
279 invertLatDes(vlistID2);
280 else
281 invertLonDes(vlistID2);
282 }
283
284 const auto streamID2 = cdo_open_write(1);
285
286 cdo_def_vlist(streamID2, vlistID2);
287
288 const auto gridsizemax = vlistGridsizeMax(vlistID1);
289 Varray<double> array1(gridsizemax), array2(gridsizemax);
290
291 int tsID = 0;
292 while (true)
293 {
294 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
295 if (nrecs == 0) break;
296
297 cdo_taxis_copy_timestep(taxisID2, taxisID1);
298
299 cdo_def_timestep(streamID2, tsID);
300
301 for (int recID = 0; recID < nrecs; recID++)
302 {
303 int varID, levelID;
304 cdo_inq_record(streamID1, &varID, &levelID);
305 size_t nmiss;
306 cdo_read_record(streamID1, array1.data(), &nmiss);
307
308 cdo_def_record(streamID2, varID, levelID);
309
310 if (operfunc1 == func_all || operfunc1 == func_fld)
311 {
312 const auto gridID1 = vlistInqVarGrid(vlistID1, varID);
313
314 if (operfunc2 == func_lat)
315 invertLatData(array1.data(), array2.data(), gridID1);
316 else
317 invertLonData(array1.data(), array2.data(), gridID1);
318
319 cdo_write_record(streamID2, array2.data(), nmiss);
320 }
321 else
322 {
323 cdo_write_record(streamID2, array1.data(), nmiss);
324 }
325 }
326 tsID++;
327 }
328
329 cdo_stream_close(streamID2);
330 cdo_stream_close(streamID1);
331
332 cdo_finish();
333
334 return nullptr;
335 }
336