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 Mrotuv mrotuv Forward rotation for MPIOM data
12 */
13
14 #include <cdi.h>
15
16 #include <mpim_grid.h>
17 #include "varray.h"
18 #include "cdo_options.h"
19 #include "process_int.h"
20 #include "cdi_lockedIO.h"
21 #include "matrix_view.h"
22
23 void
rotate_uv(Varray<double> & u_i_v,Varray<double> & v_j_v,long nx,long ny,Varray<double> & lon_v,Varray<double> & lat_v,Varray<double> & u_lon_v,Varray<double> & v_lat_v)24 rotate_uv(Varray<double> &u_i_v, Varray<double> &v_j_v, long nx, long ny, Varray<double> &lon_v, Varray<double> &lat_v,
25 Varray<double> &u_lon_v, Varray<double> &v_lat_v)
26 {
27 /*
28 in :: u_i[ny][nx], v_j[ny][nx] vector components in i-j-direction
29 in :: lat[ny][nx], lon[ny][nx] latitudes and longitudes
30 out :: u_lon[ny][nx], v_lat[ny][nx] vector components in lon-lat direction
31 */
32 constexpr double pi = 3.14159265359;
33 MatrixView<double> lon(lon_v.data(), ny, nx);
34 MatrixView<double> lat(lat_v.data(), ny, nx);
35 MatrixView<double> u_i(u_i_v.data(), ny, nx);
36 MatrixView<double> v_j(v_j_v.data(), ny, nx);
37 MatrixView<double> u_lon(u_lon_v.data(), ny, nx);
38 MatrixView<double> v_lat(v_lat_v.data(), ny, nx);
39
40 // specification whether change in sign is needed for the input arrays
41 constexpr auto change_sign_u = false;
42 constexpr auto change_sign_v = true;
43
44 // initialization
45 v_lat_v.assign(nx * ny, 0.0);
46 u_lon_v.assign(nx * ny, 0.0);
47
48 // rotation
49 for (long j = 0; j < ny; j++)
50 for (long i = 0; i < nx; i++)
51 {
52 auto ip1 = i + 1;
53 auto im1 = i - 1;
54 auto jp1 = j + 1;
55 auto jm1 = j - 1;
56 if (ip1 >= nx) ip1 = 0; // the 0-meridian
57 if (im1 < 0) im1 = nx - 1;
58 if (jp1 >= ny) jp1 = j; // treatment of the last..
59 if (jm1 < 0) jm1 = j; // .. and the fist grid-row
60
61 // difference in latitudes
62 auto dlat_i = lat[j][ip1] - lat[j][im1];
63 auto dlat_j = lat[jp1][i] - lat[jm1][i];
64
65 // difference in longitudes
66 auto dlon_i = lon[j][ip1] - lon[j][im1];
67 if (dlon_i > pi) dlon_i -= 2 * pi;
68 if (dlon_i < (-pi)) dlon_i += 2 * pi;
69 auto dlon_j = lon[jp1][i] - lon[jm1][i];
70 if (dlon_j > pi) dlon_j -= 2 * pi;
71 if (dlon_j < (-pi)) dlon_j += 2 * pi;
72
73 const auto lat_factor = std::cos(lat[j][i]);
74 dlon_i = dlon_i * lat_factor;
75 dlon_j = dlon_j * lat_factor;
76
77 // projection by scalar product
78 u_lon[j][i] = u_i[j][i] * dlon_i + v_j[j][i] * dlat_i;
79 v_lat[j][i] = u_i[j][i] * dlon_j + v_j[j][i] * dlat_j;
80
81 const auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
82 const auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
83
84 if (std::fabs(dist_i) > 0.0 && std::fabs(dist_j) > 0.0)
85 {
86 u_lon[j][i] /= dist_i;
87 v_lat[j][i] /= dist_j;
88 }
89 else
90 {
91 u_lon[j][i] = 0.0;
92 v_lat[j][i] = 0.0;
93 }
94
95 // velocity vector lengths
96 auto absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
97 auto absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
98
99 u_lon[j][i] *= absold;
100 v_lat[j][i] *= absold;
101
102 if (absnew > 0.0)
103 {
104 u_lon[j][i] /= absnew;
105 v_lat[j][i] /= absnew;
106 }
107 else
108 {
109 u_lon[j][i] = 0.0;
110 v_lat[j][i] = 0.0;
111 }
112
113 // change sign
114 if (change_sign_u) u_lon[j][i] *= -1;
115 if (change_sign_v) v_lat[j][i] *= -1;
116
117 if (Options::cdoVerbose)
118 {
119 absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
120 absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
121
122 if (i % 20 == 0 && j % 20 == 0 && absold > 0.0)
123 {
124 printf("(absold,absnew) %ld %ld %g %g %g %g %g %g\n", j + 1, i + 1, absold, absnew, u_i[j][i], v_j[j][i],
125 u_lon[j][i], v_lat[j][i]);
126
127 // test orthogonality
128 if ((dlon_i * dlon_j + dlat_j * dlat_i) > 0.1)
129 fprintf(stderr, "orthogonal? %ld %ld %g\n", j + 1, i + 1, (dlon_i * dlon_j + dlat_j * dlat_i));
130 }
131 }
132 }
133 }
134
135 void
p_to_uv_grid(long nlon,long nlat,Varray<double> & grid1x_v,Varray<double> & grid1y_v,Varray<double> & gridux_v,Varray<double> & griduy_v,Varray<double> & gridvx_v,Varray<double> & gridvy_v)136 p_to_uv_grid(long nlon, long nlat, Varray<double> &grid1x_v, Varray<double> &grid1y_v, Varray<double> &gridux_v,
137 Varray<double> &griduy_v, Varray<double> &gridvx_v, Varray<double> &gridvy_v)
138 {
139 MatrixView<double> grid1x(grid1x_v.data(), nlat, nlon);
140 MatrixView<double> grid1y(grid1y_v.data(), nlat, nlon);
141 MatrixView<double> gridux(gridux_v.data(), nlat, nlon);
142 MatrixView<double> griduy(griduy_v.data(), nlat, nlon);
143 MatrixView<double> gridvx(gridvx_v.data(), nlat, nlon);
144 MatrixView<double> gridvy(gridvy_v.data(), nlat, nlon);
145
146 // interpolate scalar to u points
147 for (long j = 0; j < nlat; j++)
148 for (long i = 0; i < nlon; i++)
149 {
150 auto ip1 = i + 1;
151 if (ip1 > nlon - 1) ip1 = 0;
152
153 gridux[j][i] = (grid1x[j][i] + grid1x[j][ip1]) * 0.5;
154 if ((grid1x[j][i] > 340 && grid1x[j][ip1] < 20) || (grid1x[j][i] < 20 && grid1x[j][ip1] > 340))
155 {
156 gridux[j][i] += (gridux[j][i] < 180) ? 180 : -180;
157 }
158
159 griduy[j][i] = (grid1y[j][i] + grid1y[j][ip1]) * 0.5;
160 }
161
162 // interpolate scalar to v points
163 for (long j = 0; j < nlat; j++)
164 for (long i = 0; i < nlon; i++)
165 {
166 auto jp1 = j + 1;
167 if (jp1 > nlat - 1) jp1 = nlat - 1;
168
169 gridvx[j][i] = (grid1x[j][i] + grid1x[jp1][i]) * 0.5;
170 if ((grid1x[j][i] > 340 && grid1x[jp1][i] < 20) || (grid1x[j][i] < 20 && grid1x[jp1][i] > 340))
171 {
172 gridvx[j][i] += (gridvx[j][i] < 180) ? 180 : -180;
173 }
174
175 gridvy[j][i] = (grid1y[j][i] + grid1y[jp1][i]) * 0.5;
176 }
177 }
178
179 void *
Mrotuv(void * process)180 Mrotuv(void *process)
181 {
182 size_t nmiss1 = 0, nmiss2 = 0;
183 int uid = -1, vid = -1;
184
185 cdo_initialize(process);
186
187 operator_check_argc(0);
188
189 const auto streamID1 = cdo_open_read(0);
190
191 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
192
193 const auto nvars = vlistNvars(vlistID1);
194 for (int varid = 0; varid < nvars; varid++)
195 {
196 const auto code = vlistInqVarCode(vlistID1, varid);
197 if (code == 3 || code == 131) uid = varid;
198 if (code == 4 || code == 132) vid = varid;
199 }
200
201 if (uid == -1 || vid == -1)
202 {
203 if (nvars == 2)
204 {
205 uid = 0;
206 vid = 1;
207 }
208 else
209 cdo_abort("U and V not found in %s", cdo_get_stream_name(0));
210 }
211
212 const auto nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, uid));
213 if (nlevs != zaxisInqSize(vlistInqVarZaxis(vlistID1, vid))) cdo_abort("U and V have different number of levels!");
214
215 auto gridID1 = vlistInqVarGrid(vlistID1, uid);
216 const auto gridID2 = vlistInqVarGrid(vlistID1, vid);
217 const auto gridsize = gridInqSize(gridID1);
218 if (gridID1 != gridID2) cdo_abort("Input grids differ!");
219
220 if (gridInqType(gridID1) != GRID_LONLAT && gridInqType(gridID1) != GRID_GAUSSIAN && gridInqType(gridID1) != GRID_CURVILINEAR)
221 cdo_abort("Grid %s unsupported!", gridNamePtr(gridInqType(gridID1)));
222
223 if (gridInqType(gridID1) != GRID_CURVILINEAR) gridID1 = gridToCurvilinear(gridID1, 0);
224
225 if (gridsize != gridInqSize(gridID1)) cdo_abort("Internal problem: gridsize changed!");
226
227 const auto nlon = gridInqXsize(gridID1);
228 const auto nlat = gridInqYsize(gridID1);
229
230 Varray<double> grid1x(gridsize), grid1y(gridsize);
231 Varray<double> gridux(gridsize), griduy(gridsize);
232 Varray<double> gridvx(gridsize), gridvy(gridsize);
233
234 gridInqXvals(gridID1, grid1x.data());
235 gridInqYvals(gridID1, grid1y.data());
236
237 // Convert lat/lon units if required
238 cdo_grid_to_degree(gridID1, CDI_XAXIS, gridsize, grid1x.data(), "grid center lon");
239 cdo_grid_to_degree(gridID1, CDI_YAXIS, gridsize, grid1y.data(), "grid center lat");
240
241 p_to_uv_grid(nlon, nlat, grid1x, grid1y, gridux, griduy, gridvx, gridvy);
242
243 const auto gridIDu = gridCreate(GRID_CURVILINEAR, nlon * nlat);
244 int datatype = CDI_UNDEFID;
245 cdiInqKeyInt(gridID1, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
246 cdiDefKeyInt(gridIDu, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
247 gridDefXsize(gridIDu, nlon);
248 gridDefYsize(gridIDu, nlat);
249 gridDefXvals(gridIDu, gridux.data());
250 gridDefYvals(gridIDu, griduy.data());
251
252 const auto gridIDv = gridCreate(GRID_CURVILINEAR, nlon * nlat);
253 cdiDefKeyInt(gridIDv, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
254 gridDefXsize(gridIDv, nlon);
255 gridDefYsize(gridIDv, nlat);
256 gridDefXvals(gridIDv, gridvx.data());
257 gridDefYvals(gridIDv, gridvy.data());
258
259 for (size_t i = 0; i < gridsize; i++)
260 {
261 grid1x[i] *= DEG2RAD;
262 grid1y[i] *= DEG2RAD;
263 }
264
265 vlistClearFlag(vlistID1);
266 for (int lid = 0; lid < nlevs; lid++) vlistDefFlag(vlistID1, uid, lid, true);
267 const auto vlistID2 = vlistCreate();
268 cdo_vlist_copy_flag(vlistID2, vlistID1);
269 vlistChangeVarGrid(vlistID2, 0, gridIDu);
270
271 vlistClearFlag(vlistID1);
272 for (int lid = 0; lid < nlevs; lid++) vlistDefFlag(vlistID1, vid, lid, true);
273 const auto vlistID3 = vlistCreate();
274 cdo_vlist_copy_flag(vlistID3, vlistID1);
275 vlistChangeVarGrid(vlistID3, 0, gridIDv);
276
277 const auto taxisID1 = vlistInqTaxis(vlistID1);
278 const auto taxisID2 = taxisDuplicate(taxisID1);
279 const auto taxisID3 = taxisDuplicate(taxisID1);
280 vlistDefTaxis(vlistID2, taxisID2);
281 vlistDefTaxis(vlistID3, taxisID3);
282
283 const auto streamID2 = cdo_open_write(1);
284 const auto streamID3 = cdo_open_write(2);
285
286 cdo_def_vlist(streamID2, vlistID2);
287 cdo_def_vlist(streamID3, vlistID3);
288
289 const auto missval1 = vlistInqVarMissval(vlistID1, uid);
290 const auto missval2 = vlistInqVarMissval(vlistID1, vid);
291
292 Varray<double> ufield_v(gridsize), vfield_v(gridsize);
293 MatrixView<double> ufield(ufield_v.data(), nlat, nlon);
294 MatrixView<double> vfield(vfield_v.data(), nlat, nlon);
295
296 Varray2D<double> urfield(nlevs), vrfield(nlevs);
297 for (int lid = 0; lid < nlevs; lid++)
298 {
299 urfield[lid].resize(gridsize);
300 vrfield[lid].resize(gridsize);
301 }
302
303 Varray2D<double> uhelp(nlat, Varray<double>(nlon + 2));
304 Varray2D<double> vhelp(nlat, Varray<double>(nlon + 2));
305
306 int tsID = 0;
307 while (true)
308 {
309 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
310 if (nrecs == 0) break;
311
312 cdo_taxis_copy_timestep(taxisID2, taxisID1);
313 cdo_def_timestep(streamID2, tsID);
314 cdo_taxis_copy_timestep(taxisID3, taxisID1);
315 cdo_def_timestep(streamID3, tsID);
316
317 for (int recID = 0; recID < nrecs; recID++)
318 {
319 int varID, levelID;
320 cdo_inq_record(streamID1, &varID, &levelID);
321
322 if (varID == uid) cdo_read_record(streamID1, urfield[levelID].data(), &nmiss1);
323 if (varID == vid) cdo_read_record(streamID1, vrfield[levelID].data(), &nmiss2);
324 }
325
326 for (int levelID = 0; levelID < nlevs; levelID++)
327 {
328 // remove missing values
329 if (nmiss1 || nmiss2)
330 {
331 for (size_t i = 0; i < gridsize; i++)
332 {
333 if (DBL_IS_EQUAL(urfield[levelID][i], missval1)) urfield[levelID][i] = 0.0;
334 if (DBL_IS_EQUAL(vrfield[levelID][i], missval2)) vrfield[levelID][i] = 0.0;
335 }
336 }
337
338 // rotate
339 rotate_uv(urfield[levelID], vrfield[levelID], nlon, nlat, grid1x, grid1y, ufield_v, vfield_v);
340
341 // load to a help field
342 for (size_t j = 0; j < nlat; j++)
343 for (size_t i = 0; i < nlon; i++)
344 {
345 uhelp[j][i + 1] = ufield[j][i];
346 vhelp[j][i + 1] = vfield[j][i];
347 }
348
349 // make help field cyclic
350 for (size_t j = 0; j < nlat; j++)
351 {
352 uhelp[j][0] = uhelp[j][nlon];
353 uhelp[j][nlon + 1] = uhelp[j][1];
354 vhelp[j][0] = vhelp[j][nlon];
355 vhelp[j][nlon + 1] = vhelp[j][1];
356 }
357
358 // interpolate on u/v points
359 for (size_t j = 0; j < nlat; j++)
360 for (size_t i = 0; i < nlon; i++)
361 {
362 ufield[j][i] = (uhelp[j][i + 1] + uhelp[j][i + 2]) * 0.5;
363 }
364
365 for (size_t j = 0; j < nlat - 1; j++)
366 for (size_t i = 0; i < nlon; i++)
367 {
368 vfield[j][i] = (vhelp[j][i + 1] + vhelp[j + 1][i + 1]) * 0.5;
369 }
370
371 for (size_t i = 0; i < nlon; i++)
372 {
373 vfield[nlat - 1][i] = vhelp[nlat - 1][i + 1];
374 }
375
376 cdo_def_record(streamID2, 0, levelID);
377 cdo_write_record(streamID2, ufield_v.data(), nmiss1);
378 cdo_def_record(streamID3, 0, levelID);
379 cdo_write_record(streamID3, vfield_v.data(), nmiss2);
380 }
381
382 tsID++;
383 }
384
385 cdo_stream_close(streamID3);
386 cdo_stream_close(streamID2);
387 cdo_stream_close(streamID1);
388
389 cdo_finish();
390
391 return nullptr;
392 }
393