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 Mrotuvb mrotuvb Backward rotation for MPIOM data
12 */
13
14 #include <cdi.h>
15
16 #include "varray.h"
17 #include "cdo_options.h"
18 #include "process_int.h"
19 #include <mpim_grid.h>
20 #include "matrix_view.h"
21
22 /*
23 !----------------------------------------------------------------------
24 !
25 ! rotation of vectors: In ocean models with rotated grids velocity
26 ! vectors are given in the direction of grid lines and rows. They
27 ! have to be rotated in latitudinal and longitudinal direction.
28 !
29 ! Note: This routine assumes positive meridional flow for a flow
30 ! from grid point(i,j) to grid point(i,j+1) and positive
31 ! zonal flow for a flow from grid point(i,j) to point(i+1,j).
32 ! this is not the case for mpi-om!
33 !
34 ! If this routine is used to rotate data of mpi-om, the
35 ! logical change_sign_v needs to be true.
36 !
37 ! j. jungclaus: 22.01.04:
38 ! Note here for the coupling fields u_i,v_j are on the non-verlapping
39 ! (ie-2=ix) grid, furthermore, the velocity fields were previously
40 ! interpolated onto the scalar points!
41 !----------------------------------------------------------------------
42 */
43
44 void
rotate_uv2(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)45 rotate_uv2(Varray<double> &u_i_v, Varray<double> &v_j_v, long nx, long ny, Varray<double> &lon_v, Varray<double> &lat_v,
46 Varray<double> &u_lon_v, Varray<double> &v_lat_v)
47 {
48 /*
49 in :: u_i[ny][nx], v_j[ny][nx] vector components in i-j-direction
50 in :: lat[ny][nx], lon[ny][nx] latitudes and longitudes
51 out :: u_lon[ny][nx], v_lat[ny][nx] vector components in lon-lat direction
52 */
53 constexpr double pi = 3.14159265359;
54 MatrixView<double> lon(lon_v.data(), ny, nx);
55 MatrixView<double> lat(lat_v.data(), ny, nx);
56 MatrixView<double> u_i(u_i_v.data(), ny, nx);
57 MatrixView<double> v_j(v_j_v.data(), ny, nx);
58 MatrixView<double> u_lon(u_lon_v.data(), ny, nx);
59 MatrixView<double> v_lat(v_lat_v.data(), ny, nx);
60
61 // specification whether change in sign is needed for the input arrays
62 constexpr auto change_sign_u = false;
63 constexpr auto change_sign_v = true;
64
65 // initialization
66 v_lat_v.assign(nx * ny, 0.0);
67 u_lon_v.assign(nx * ny, 0.0);
68
69 // change sign
70 if (change_sign_u)
71 for (long i = 0; i < nx * ny; i++) u_i_v[i] *= -1;
72
73 if (change_sign_v)
74 for (long i = 0; i < nx * ny; i++) v_j_v[i] *= -1;
75
76 // rotation
77 for (long j = 0; j < ny; j++)
78 for (long i = 0; i < nx; i++)
79 {
80 auto ip1 = i + 1;
81 auto im1 = i - 1;
82 auto jp1 = j + 1;
83 auto jm1 = j - 1;
84 if (ip1 >= nx) ip1 = 0; // the 0-meridian
85 if (im1 < 0) im1 = nx - 1;
86 if (jp1 >= ny) jp1 = j; // treatment of the last..
87 if (jm1 < 0) jm1 = j; // .. and the fist grid-row
88
89 // difference in latitudes
90 auto dlat_i = lat[j][ip1] - lat[j][im1];
91 auto dlat_j = lat[jp1][i] - lat[jm1][i];
92
93 // difference in longitudes
94 auto dlon_i = lon[j][ip1] - lon[j][im1];
95 if (dlon_i > pi) dlon_i -= 2 * pi;
96 if (dlon_i < (-pi)) dlon_i += 2 * pi;
97 auto dlon_j = lon[jp1][i] - lon[jm1][i];
98 if (dlon_j > pi) dlon_j -= 2 * pi;
99 if (dlon_j < (-pi)) dlon_j += 2 * pi;
100
101 const auto lat_factor = std::cos(lat[j][i]);
102 dlon_i = dlon_i * lat_factor;
103 dlon_j = dlon_j * lat_factor;
104
105 // projection by scalar product
106 u_lon[j][i] = u_i[j][i] * dlon_i + v_j[j][i] * dlat_i;
107 v_lat[j][i] = u_i[j][i] * dlon_j + v_j[j][i] * dlat_j;
108
109 const auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
110 const auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
111
112 if (std::fabs(dist_i) > 0 && std::fabs(dist_j) > 0)
113 {
114 u_lon[j][i] /= dist_i;
115 v_lat[j][i] /= dist_j;
116 }
117 else
118 {
119 u_lon[j][i] = 0.0;
120 v_lat[j][i] = 0.0;
121 }
122
123 if (Options::cdoVerbose)
124 {
125 // velocity vector lengths
126 const auto absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
127 const auto absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
128
129 if (i % 20 == 0 && j % 20 == 0 && absold > 0)
130 {
131 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],
132 u_lon[j][i], v_lat[j][i]);
133
134 // test orthogonality
135 if ((dlon_i * dlon_j + dlat_j * dlat_i) > 0.1)
136 fprintf(stderr, "orthogonal? %ld %ld %g\n", j + 1, i + 1, (dlon_i * dlon_j + dlat_j * dlat_i));
137 }
138 }
139 }
140 }
141
142 static void
uv_to_p_grid(size_t nlon,size_t nlat,Varray<double> & grid1x_v,Varray<double> & grid1y_v,Varray<double> & grid2x_v,Varray<double> & grid2y_v,Varray<double> & grid3x_v,Varray<double> & grid3y_v)143 uv_to_p_grid(size_t nlon, size_t nlat, Varray<double> &grid1x_v, Varray<double> &grid1y_v, Varray<double> &grid2x_v,
144 Varray<double> &grid2y_v, Varray<double> &grid3x_v, Varray<double> &grid3y_v)
145 {
146 MatrixView<double> grid1x(grid1x_v.data(), nlat, nlon);
147 MatrixView<double> grid1y(grid1y_v.data(), nlat, nlon);
148 MatrixView<double> grid2x(grid2x_v.data(), nlat, nlon);
149 MatrixView<double> grid2y(grid2y_v.data(), nlat, nlon);
150 MatrixView<double> grid3x(grid3x_v.data(), nlat, nlon);
151 MatrixView<double> grid3y(grid3y_v.data(), nlat, nlon);
152
153 Varray2D<double> gxhelp(nlat, Varray<double>(nlon + 2)), gyhelp(nlat, Varray<double>(nlon + 2));
154
155 // load to a help field
156 for (size_t j = 0; j < nlat; j++)
157 for (size_t i = 0; i < nlon; i++)
158 {
159 gxhelp[j][i + 1] = grid1x[j][i];
160 gyhelp[j][i + 1] = grid1y[j][i];
161 }
162
163 // make help field cyclic
164 for (size_t j = 0; j < nlat; j++)
165 {
166 gxhelp[j][0] = gxhelp[j][nlon];
167 gxhelp[j][nlon + 1] = gxhelp[j][1];
168 gyhelp[j][0] = gyhelp[j][nlon];
169 gyhelp[j][nlon + 1] = gyhelp[j][1];
170 }
171
172 // interpolate u to scalar points
173 for (size_t j = 0; j < nlat; j++)
174 for (size_t i = 0; i < nlon; i++)
175 {
176 grid3x[j][i] = (gxhelp[j][i] + gxhelp[j][i + 1]) * 0.5;
177 if ((gxhelp[j][i] > 340 && gxhelp[j][i + 1] < 20) || (gxhelp[j][i] < 20 && gxhelp[j][i + 1] > 340))
178 {
179 grid3x[j][i] += (grid3x[j][i] < 180) ? 180 : -180;
180 }
181
182 grid3y[j][i] = (gyhelp[j][i] + gyhelp[j][i + 1]) * 0.5;
183 }
184
185 // load to a help field
186 for (size_t j = 0; j < nlat; j++)
187 for (size_t i = 0; i < nlon; i++)
188 {
189 gxhelp[j][i + 1] = grid2x[j][i];
190 gyhelp[j][i + 1] = grid2y[j][i];
191 }
192
193 // make help field cyclic
194 for (size_t j = 0; j < nlat; j++)
195 {
196 gxhelp[j][0] = gxhelp[j][nlon];
197 gxhelp[j][nlon + 1] = gxhelp[j][1];
198 gyhelp[j][0] = gyhelp[j][nlon];
199 gyhelp[j][nlon + 1] = gyhelp[j][1];
200 }
201
202 // interpolate v to scalar points
203 for (size_t j = 1; j < nlat - 1; j++)
204 for (size_t i = 0; i < nlon; i++)
205 {
206 auto gx = (gxhelp[j][i + 1] + gxhelp[j - 1][i + 1]) * 0.5;
207 if ((gxhelp[j][i + 1] > 340 && gxhelp[j - 1][i + 1] < 20) || (gxhelp[j][i + 1] < 20 && gxhelp[j - 1][i + 1] > 340))
208 {
209 gx += (gx < 180) ? 180 : -180;
210 }
211
212 auto gy = (gyhelp[j][i + 1] + gyhelp[j - 1][i + 1]) * 0.5;
213
214 // printf("%d %d %g %g %g %g \n", j, i, gx, gy, grid3x[j][i], grid3y[j][i]);
215
216 auto gx2 = (gx + grid3x[j][i]) * 0.5;
217 if ((gx > 340 && grid3x[j][i] < 20) || (gx < 20 && grid3x[j][i] > 340))
218 {
219 gx2 += (gx2 < 180) ? 180 : -180;
220 }
221
222 auto gy2 = (gy + grid3y[j][i]) * 0.5;
223
224 grid3x[j][i] = gx2;
225 grid3y[j][i] = gy2;
226
227 // printf("%d %d %g %g %g %g \n", j, i, gx2, gy2, grid3x[j][i], grid3y[j][i]);
228 }
229 }
230
231 void *
Mrotuvb(void * process)232 Mrotuvb(void *process)
233 {
234 cdo_initialize(process);
235
236 const auto gpint = !(cdo_operator_argc() == 1 && cdo_operator_argv(0) == "noint");
237
238 const auto streamID1 = cdo_open_read(0);
239 const auto streamID2 = cdo_open_read(1);
240
241 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
242 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
243
244 auto nvars = vlistNvars(vlistID1);
245 if (nvars > 1) cdo_abort("More than one variable found in %s", cdo_get_stream_name(0));
246 nvars = vlistNvars(vlistID2);
247 if (nvars > 1) cdo_abort("More than one variable found in %s", cdo_get_stream_name(1));
248
249 auto gridID1 = vlistGrid(vlistID1, 0);
250 auto gridID2 = vlistGrid(vlistID2, 0);
251 const auto gridsize = gridInqSize(gridID1);
252 if (gpint && gridID1 == gridID2) cdo_abort("Input grids are the same, use parameter >noint< to disable interpolation!");
253 if (!gpint && gridID1 != gridID2) cdo_abort("Input grids are not the same!");
254 if (gridsize != gridInqSize(gridID2)) cdo_abort("Grids have different size!");
255
256 if (gridInqType(gridID1) != GRID_LONLAT && gridInqType(gridID1) != GRID_GAUSSIAN && gridInqType(gridID1) != GRID_CURVILINEAR)
257 cdo_abort("Grid %s unsupported!", gridNamePtr(gridInqType(gridID1)));
258
259 if (gridInqType(gridID1) != GRID_CURVILINEAR) gridID1 = gridToCurvilinear(gridID1, 1);
260
261 if (gridsize != gridInqSize(gridID1)) cdo_abort("Internal problem: gridsize changed!");
262
263 if (gridInqType(gridID2) != GRID_CURVILINEAR) gridID2 = gridToCurvilinear(gridID2, 1);
264
265 if (gridsize != gridInqSize(gridID2)) cdo_abort("Internal problem: gridsize changed!");
266
267 const auto nlon = gridInqXsize(gridID1);
268 const auto nlat = gridInqYsize(gridID1);
269
270 Varray<double> grid1x(gridsize), grid1y(gridsize);
271 Varray<double> grid2x(gridsize), grid2y(gridsize);
272 Varray<double> grid3x(gridsize), grid3y(gridsize);
273
274 gridInqXvals(gridID1, grid1x.data());
275 gridInqYvals(gridID1, grid1y.data());
276
277 // Convert lat/lon units if required
278 cdo_grid_to_degree(gridID1, CDI_XAXIS, gridsize, grid1x.data(), "grid1 center lon");
279 cdo_grid_to_degree(gridID1, CDI_YAXIS, gridsize, grid1y.data(), "grid1 center lat");
280
281 gridInqXvals(gridID2, grid2x.data());
282 gridInqYvals(gridID2, grid2y.data());
283
284 // Convert lat/lon units if required
285 cdo_grid_to_degree(gridID2, CDI_XAXIS, gridsize, grid2x.data(), "grid2 center lon");
286 cdo_grid_to_degree(gridID2, CDI_YAXIS, gridsize, grid2y.data(), "grid2 center lat");
287
288 if (gpint)
289 {
290 uv_to_p_grid(nlon, nlat, grid1x, grid1y, grid2x, grid2y, grid3x, grid3y);
291 }
292 else
293 {
294 grid3x = grid1x;
295 grid3y = grid1y;
296 }
297
298 const auto gridID3 = gridCreate(GRID_CURVILINEAR, gridsize);
299 int datatype = CDI_UNDEFID;
300 cdiInqKeyInt(gridID1, CDI_GLOBAL, CDI_KEY_DATATYPE, &datatype);
301 cdiDefKeyInt(gridID3, CDI_GLOBAL, CDI_KEY_DATATYPE, datatype);
302 gridDefXsize(gridID3, nlon);
303 gridDefYsize(gridID3, nlat);
304 gridDefXvals(gridID3, grid3x.data());
305 gridDefYvals(gridID3, grid3y.data());
306
307 for (size_t i = 0; i < gridsize; i++)
308 {
309 grid3x[i] *= DEG2RAD;
310 grid3y[i] *= DEG2RAD;
311 }
312
313 const auto vlistID3 = vlistCreate();
314 vlistCopy(vlistID3, vlistID1);
315 vlistCat(vlistID3, vlistID2);
316
317 const auto code1 = vlistInqVarCode(vlistID1, 0);
318 const auto code2 = vlistInqVarCode(vlistID2, 0);
319
320 if (code1 == code2) vlistDefVarCode(vlistID3, 1, code1 + 1);
321
322 vlistChangeGrid(vlistID3, gridID1, gridID3);
323 vlistChangeGrid(vlistID3, gridID2, gridID3);
324
325 const auto taxisID1 = vlistInqTaxis(vlistID1);
326 const auto taxisID3 = taxisDuplicate(taxisID1);
327 vlistDefTaxis(vlistID3, taxisID3);
328
329 if (Options::cdoVerbose) vlistPrint(vlistID3);
330
331 const auto streamID3 = cdo_open_write(2);
332 cdo_def_vlist(streamID3, vlistID3);
333
334 const auto missval1 = vlistInqVarMissval(vlistID1, 0);
335 const auto missval2 = vlistInqVarMissval(vlistID2, 0);
336
337 Varray<double> urfield(gridsize), vrfield(gridsize);
338 Varray<double> ufield_v(gridsize), vfield_v(gridsize);
339 MatrixView<double> ufield(ufield_v.data(), nlat, nlon);
340 MatrixView<double> vfield(vfield_v.data(), nlat, nlon);
341
342 Varray<double> uhelp_v, vhelp_v;
343 if (gpint)
344 {
345 const auto gridsizex = (nlon + 2) * nlat;
346 uhelp_v.resize(gridsizex);
347 vhelp_v.resize(gridsizex);
348 }
349 MatrixView<double> uhelp(uhelp_v.data(), nlat, nlon + 2);
350 MatrixView<double> vhelp(vhelp_v.data(), nlat, nlon + 2);
351
352 int tsID = 0;
353 while (true)
354 {
355 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
356 if (nrecs == 0) break;
357
358 cdo_taxis_copy_timestep(taxisID3, taxisID1);
359
360 cdo_def_timestep(streamID3, tsID);
361
362 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
363
364 if (nrecs != nrecs2) cdo_abort("Input streams have different number of levels!");
365
366 for (int recID = 0; recID < nrecs; recID++)
367 {
368 int varID1, varID2, levelID;
369 cdo_inq_record(streamID1, &varID1, &levelID);
370 cdo_inq_record(streamID2, &varID2, &levelID);
371
372 size_t nmiss1, nmiss2;
373 cdo_read_record(streamID1, ufield_v.data(), &nmiss1);
374 cdo_read_record(streamID2, vfield_v.data(), &nmiss2);
375
376 // remove missing values
377 if (nmiss1 || nmiss2)
378 {
379 for (size_t i = 0; i < gridsize; i++)
380 {
381 if (DBL_IS_EQUAL(ufield_v[i], missval1)) ufield_v[i] = 0.0;
382 if (DBL_IS_EQUAL(vfield_v[i], missval2)) vfield_v[i] = 0.0;
383 }
384 }
385
386 if (gpint)
387 {
388 // load to a help field
389 for (size_t j = 0; j < nlat; j++)
390 for (size_t i = 0; i < nlon; i++)
391 {
392 uhelp[j][i + 1] = ufield[j][i];
393 vhelp[j][i + 1] = vfield[j][i];
394 }
395
396 // make help field cyclic
397 for (size_t j = 0; j < nlat; j++)
398 {
399 uhelp[j][0] = uhelp[j][nlon];
400 uhelp[j][nlon + 1] = uhelp[j][1];
401 vhelp[j][0] = vhelp[j][nlon];
402 vhelp[j][nlon + 1] = vhelp[j][1];
403 }
404
405 // interpolate on pressure points
406 for (size_t j = 1; j < nlat; j++)
407 for (size_t i = 0; i < nlon; i++)
408 {
409 ufield[j][i] = (uhelp[j][i] + uhelp[j][i + 1]) * 0.5;
410 vfield[j][i] = (vhelp[j - 1][i + 1] + vhelp[j][i + 1]) * 0.5;
411 }
412 }
413
414 for (size_t i = 0; i < nlon; i++)
415 {
416 ufield[0][i] = 0.0;
417 vfield[0][i] = 0.0;
418 }
419
420 // rotate
421 rotate_uv2(ufield_v, vfield_v, nlon, nlat, grid3x, grid3y, urfield, vrfield);
422
423 // calc lat, lon, Auv and alpha
424 /*
425 {
426 double lat, lon, auv, alpha;
427 for ( int j = 1; j < nlat-1; j += 3 )
428 for ( int i = 0; i < nlon; i += 3 )
429 {
430 lat = grid3y[j][i]*RAD2DEG;
431 lon = grid3x[j][i]*RAD2DEG;
432 auv = std::sqrt(urfield[j][i]*urfield[j][i] + vrfield[j][i]*vrfield[j][i]);
433 alpha = std::atan2(vrfield[j][i], urfield[j][i]);
434 alpha = 90. - alpha*RAD2DEG;
435
436 if ( alpha < 0 ) alpha += 360.;
437 if ( alpha > 360 ) alpha -= 360.;
438
439 printf("%g %g %g %g\n", lon, lat, alpha, auv);
440 }
441 }
442 */
443 cdo_def_record(streamID3, 0, levelID);
444 cdo_write_record(streamID3, urfield.data(), 0);
445 cdo_def_record(streamID3, 1, levelID);
446 cdo_write_record(streamID3, vrfield.data(), 0);
447 }
448
449 tsID++;
450 }
451
452 cdo_stream_close(streamID3);
453 cdo_stream_close(streamID2);
454 cdo_stream_close(streamID1);
455
456 cdo_finish();
457
458 return nullptr;
459 }
460