Commit cc4b668a authored by Uwe Schulzweida's avatar Uwe Schulzweida

Replaced IX2D() by matrix_view.

parent 9f5e8235
Pipeline #6840 failed with stages
in 5 minutes and 51 seconds
......@@ -231,8 +231,8 @@ fillmiss_one_step(Field &field1, Field &field2, int maxfill)
kv = ko + ku;
if (kh == 0)
{
s1 = 0.;
k1 = 0;
s1 = 0.0;
k1 = 0.0;
}
else if (kl == 0)
{
......@@ -295,20 +295,11 @@ fillmiss_one_step(Field &field1, Field &field2, int maxfill)
else if (k2 == 0)
matrix2[j][i] = s1;
else
{
if (k1 <= k2)
{
matrix2[j][i] = s1;
}
else
{
matrix2[j][i] = s2;
}
}
matrix2[j][i] = (k1 <= k2) ? s1 : s2;
// printf("%d %d %2d %2d %2d %2d %2g %2g %2g %2g %2g %2g %2g\n",
// j,i,kr,kl,ku,ko,xr,xl,xu,xo,s1,s2,matrix2[j][i]);
/* matrix1[j][i] = matrix2[j][i]; */
// matrix1[j][i] = matrix2[j][i];
}
else
{
......@@ -329,7 +320,6 @@ setmisstodis(Field &field1, Field &field2, int numNeighbors)
auto array1 = field1.vec_d.data();
auto array2 = field2.vec_d.data();
const auto gridtype = gridInqType(gridID);
const auto gridsize = gridInqSize(gridID);
auto nmiss = field1.nmiss;
......
......@@ -28,106 +28,111 @@
#include "cdo_options.h"
#include "process_int.h"
#include "cdi_lockedIO.h"
#include "matrix_view.h"
void
rotate_uv(double *u_i, double *v_j, long ix, long iy, double *lon, double *lat, double *u_lon, double *v_lat)
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)
{
/*
in :: u_i(ix,iy),v_j(ix,iy) vector components in i-j-direction
in :: lat(ix,iy),lon(ix,iy) latitudes and longitudes
out :: u_lon(ix,iy),v_lat(ix,iy) vector components in lon-lat direction
in :: u_i[ny][nx], v_j[ny][nx] vector components in i-j-direction
in :: lat[ny][nx], lon[ny][nx] latitudes and longitudes
out :: u_lon[ny][nx], v_lat[ny][nx] vector components in lon-lat direction
*/
double pi = 3.14159265359;
constexpr double pi = 3.14159265359;
MatrixView<double> lon(lon_v.data(), ny, nx);
MatrixView<double> lat(lat_v.data(), ny, nx);
MatrixView<double> u_i(u_i_v.data(), ny, nx);
MatrixView<double> v_j(v_j_v.data(), ny, nx);
MatrixView<double> u_lon(u_lon_v.data(), ny, nx);
MatrixView<double> v_lat(v_lat_v.data(), ny, nx);
// specification whether change in sign is needed for the input arrays
bool change_sign_u = false;
bool change_sign_v = true;
constexpr bool change_sign_u = false;
constexpr bool change_sign_v = true;
// initialization
for (long i = 0; i < ix * iy; i++)
{
v_lat[i] = 0;
u_lon[i] = 0;
}
v_lat_v.assign(nx * ny, 0.0);
u_lon_v.assign(nx * ny, 0.0);
// rotation
for (long j = 0; j < iy; j++)
for (long i = 0; i < ix; i++)
for (long j = 0; j < ny; j++)
for (long i = 0; i < nx; i++)
{
auto ip1 = i + 1;
auto im1 = i - 1;
auto jp1 = j + 1;
auto jm1 = j - 1;
if (ip1 >= ix) ip1 = 0; // the 0-meridian
if (im1 < 0) im1 = ix - 1;
if (jp1 >= iy) jp1 = j; // treatment of the last..
if (ip1 >= nx) ip1 = 0; // the 0-meridian
if (im1 < 0) im1 = nx - 1;
if (jp1 >= ny) jp1 = j; // treatment of the last..
if (jm1 < 0) jm1 = j; // .. and the fist grid-row
// difference in latitudes
auto dlat_i = lat[IX2D(j, ip1, ix)] - lat[IX2D(j, im1, ix)];
auto dlat_j = lat[IX2D(jp1, i, ix)] - lat[IX2D(jm1, i, ix)];
auto dlat_i = lat[j][ip1] - lat[j][im1];
auto dlat_j = lat[jp1][i] - lat[jm1][i];
// difference in longitudes
auto dlon_i = lon[IX2D(j, ip1, ix)] - lon[IX2D(j, im1, ix)];
auto dlon_i = lon[j][ip1] - lon[j][im1];
if (dlon_i > pi) dlon_i -= 2 * pi;
if (dlon_i < (-pi)) dlon_i += 2 * pi;
auto dlon_j = lon[IX2D(jp1, i, ix)] - lon[IX2D(jm1, i, ix)];
auto dlon_j = lon[jp1][i] - lon[jm1][i];
if (dlon_j > pi) dlon_j -= 2 * pi;
if (dlon_j < (-pi)) dlon_j += 2 * pi;
const auto lat_factor = std::cos(lat[IX2D(j, i, ix)]);
const auto lat_factor = std::cos(lat[j][i]);
dlon_i = dlon_i * lat_factor;
dlon_j = dlon_j * lat_factor;
// projection by scalar product
u_lon[IX2D(j, i, ix)] = u_i[IX2D(j, i, ix)] * dlon_i + v_j[IX2D(j, i, ix)] * dlat_i;
v_lat[IX2D(j, i, ix)] = u_i[IX2D(j, i, ix)] * dlon_j + v_j[IX2D(j, i, ix)] * dlat_j;
u_lon[j][i] = u_i[j][i] * dlon_i + v_j[j][i] * dlat_i;
v_lat[j][i] = u_i[j][i] * dlon_j + v_j[j][i] * dlat_j;
auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
const auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
const auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
if (std::fabs(dist_i) > 0 && std::fabs(dist_j) > 0)
if (std::fabs(dist_i) > 0.0 && std::fabs(dist_j) > 0.0)
{
u_lon[IX2D(j, i, ix)] /= dist_i;
v_lat[IX2D(j, i, ix)] /= dist_j;
u_lon[j][i] /= dist_i;
v_lat[j][i] /= dist_j;
}
else
{
u_lon[IX2D(j, i, ix)] = 0.0;
v_lat[IX2D(j, i, ix)] = 0.0;
u_lon[j][i] = 0.0;
v_lat[j][i] = 0.0;
}
// velocity vector lengths
auto absold = std::sqrt(u_i[IX2D(j, i, ix)] * u_i[IX2D(j, i, ix)] + v_j[IX2D(j, i, ix)] * v_j[IX2D(j, i, ix)]);
auto absnew = std::sqrt(u_lon[IX2D(j, i, ix)] * u_lon[IX2D(j, i, ix)] + v_lat[IX2D(j, i, ix)] * v_lat[IX2D(j, i, ix)]);
auto absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
auto absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
u_lon[IX2D(j, i, ix)] *= absold;
v_lat[IX2D(j, i, ix)] *= absold;
u_lon[j][i] *= absold;
v_lat[j][i] *= absold;
if (absnew > 0)
if (absnew > 0.0)
{
u_lon[IX2D(j, i, ix)] /= absnew;
v_lat[IX2D(j, i, ix)] /= absnew;
u_lon[j][i] /= absnew;
v_lat[j][i] /= absnew;
}
else
{
u_lon[IX2D(j, i, ix)] = 0.0;
v_lat[IX2D(j, i, ix)] = 0.0;
u_lon[j][i] = 0.0;
v_lat[j][i] = 0.0;
}
// change sign
if (change_sign_u) u_lon[IX2D(j, i, ix)] *= -1;
if (change_sign_v) v_lat[IX2D(j, i, ix)] *= -1;
if (change_sign_u) u_lon[j][i] *= -1;
if (change_sign_v) v_lat[j][i] *= -1;
if (Options::cdoVerbose)
{
absold = std::sqrt(u_i[IX2D(j, i, ix)] * u_i[IX2D(j, i, ix)] + v_j[IX2D(j, i, ix)] * v_j[IX2D(j, i, ix)]);
absnew = std::sqrt(u_lon[IX2D(j, i, ix)] * u_lon[IX2D(j, i, ix)] + v_lat[IX2D(j, i, ix)] * v_lat[IX2D(j, i, ix)]);
absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
if (i % 20 == 0 && j % 20 == 0 && absold > 0)
if (i % 20 == 0 && j % 20 == 0 && absold > 0.0)
{
printf("(absold,absnew) %ld %ld %g %g %g %g %g %g\n", j + 1, i + 1, absold, absnew, u_i[IX2D(j, i, ix)],
v_j[IX2D(j, i, ix)], u_lon[IX2D(j, i, ix)], v_lat[IX2D(j, i, ix)]);
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], u_lon[j][i], v_lat[j][i]);
// test orthogonality
if ((dlon_i * dlon_j + dlat_j * dlat_i) > 0.1)
......@@ -138,8 +143,16 @@ rotate_uv(double *u_i, double *v_j, long ix, long iy, double *lon, double *lat,
}
void
p_to_uv_grid(long nlon, long nlat, double *grid1x, double *grid1y, double *gridux, double *griduy, double *gridvx, double *gridvy)
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)
{
MatrixView<double> grid1x(grid1x_v.data(), nlat, nlon);
MatrixView<double> grid1y(grid1y_v.data(), nlat, nlon);
MatrixView<double> gridux(gridux_v.data(), nlat, nlon);
MatrixView<double> griduy(griduy_v.data(), nlat, nlon);
MatrixView<double> gridvx(gridvx_v.data(), nlat, nlon);
MatrixView<double> gridvy(gridvy_v.data(), nlat, nlon);
// interpolate scalar to u points
for (long j = 0; j < nlat; j++)
for (long i = 0; i < nlon; i++)
......@@ -147,17 +160,13 @@ p_to_uv_grid(long nlon, long nlat, double *grid1x, double *grid1y, double *gridu
auto ip1 = i + 1;
if (ip1 > nlon - 1) ip1 = 0;
gridux[IX2D(j, i, nlon)] = (grid1x[IX2D(j, i, nlon)] + grid1x[IX2D(j, ip1, nlon)]) * 0.5;
if ((grid1x[IX2D(j, i, nlon)] > 340 && grid1x[IX2D(j, ip1, nlon)] < 20)
|| (grid1x[IX2D(j, i, nlon)] < 20 && grid1x[IX2D(j, ip1, nlon)] > 340))
gridux[j][i] = (grid1x[j][i] + grid1x[j][ip1]) * 0.5;
if ((grid1x[j][i] > 340 && grid1x[j][ip1] < 20) || (grid1x[j][i] < 20 && grid1x[j][ip1] > 340))
{
if (gridux[IX2D(j, i, nlon)] < 180)
gridux[IX2D(j, i, nlon)] += 180;
else
gridux[IX2D(j, i, nlon)] -= 180;
gridux[j][i] += (gridux[j][i] < 180) ? 180 : -180;
}
griduy[IX2D(j, i, nlon)] = (grid1y[IX2D(j, i, nlon)] + grid1y[IX2D(j, ip1, nlon)]) * 0.5;
griduy[j][i] = (grid1y[j][i] + grid1y[j][ip1]) * 0.5;
}
// interpolate scalar to v points
......@@ -167,26 +176,19 @@ p_to_uv_grid(long nlon, long nlat, double *grid1x, double *grid1y, double *gridu
auto jp1 = j + 1;
if (jp1 > nlat - 1) jp1 = nlat - 1;
gridvx[IX2D(j, i, nlon)] = (grid1x[IX2D(j, i, nlon)] + grid1x[IX2D(jp1, i, nlon)]) * 0.5;
if ((grid1x[IX2D(j, i, nlon)] > 340 && grid1x[IX2D(jp1, i, nlon)] < 20)
|| (grid1x[IX2D(j, i, nlon)] < 20 && grid1x[IX2D(jp1, i, nlon)] > 340))
gridvx[j][i] = (grid1x[j][i] + grid1x[jp1][i]) * 0.5;
if ((grid1x[j][i] > 340 && grid1x[jp1][i] < 20) || (grid1x[j][i] < 20 && grid1x[jp1][i] > 340))
{
if (gridvx[IX2D(j, i, nlon)] < 180)
gridvx[IX2D(j, i, nlon)] += 180;
else
gridvx[IX2D(j, i, nlon)] -= 180;
gridvx[j][i] += (gridvx[j][i] < 180) ? 180 : -180;
}
gridvy[IX2D(j, i, nlon)] = (grid1y[IX2D(j, i, nlon)] + grid1y[IX2D(jp1, i, nlon)]) * 0.5;
gridvy[j][i] = (grid1y[j][i] + grid1y[jp1][i]) * 0.5;
}
}
void *
Mrotuv(void *process)
{
int nrecs;
int levelID;
int varID, varid;
size_t nmiss1 = 0, nmiss2 = 0;
int uid = -1, vid = -1;
......@@ -199,9 +201,9 @@ Mrotuv(void *process)
const auto vlistID1 = cdoStreamInqVlist(streamID1);
const auto nvars = vlistNvars(vlistID1);
for (varid = 0; varid < nvars; varid++)
for (int varid = 0; varid < nvars; varid++)
{
int code = vlistInqVarCode(vlistID1, varid);
const auto code = vlistInqVarCode(vlistID1, varid);
if (code == 3 || code == 131) uid = varid;
if (code == 4 || code == 132) vid = varid;
}
......@@ -239,8 +241,6 @@ Mrotuv(void *process)
Varray<double> gridux(gridsize), griduy(gridsize);
Varray<double> gridvx(gridsize), gridvy(gridsize);
const auto gridsizex = (nlon + 2) * nlat;
gridInqXvals(gridID1, grid1x.data());
gridInqYvals(gridID1, grid1y.data());
......@@ -248,7 +248,7 @@ Mrotuv(void *process)
cdo_grid_to_degree(gridID1, CDI_XAXIS, gridsize, grid1x.data(), "grid center lon");
cdo_grid_to_degree(gridID1, CDI_YAXIS, gridsize, grid1y.data(), "grid center lat");
p_to_uv_grid(nlon, nlat, grid1x.data(), grid1y.data(), gridux.data(), griduy.data(), gridvx.data(), gridvy.data());
p_to_uv_grid(nlon, nlat, grid1x, grid1y, gridux, griduy, gridvx, gridvy);
const auto gridIDu = gridCreate(GRID_CURVILINEAR, nlon * nlat);
gridDefDatatype(gridIDu, gridInqDatatype(gridID1));
......@@ -297,7 +297,9 @@ Mrotuv(void *process)
const auto missval1 = vlistInqVarMissval(vlistID1, uid);
const auto missval2 = vlistInqVarMissval(vlistID1, vid);
Varray<double> ufield(gridsize), vfield(gridsize);
Varray<double> ufield_v(gridsize), vfield_v(gridsize);
MatrixView<double> ufield(ufield_v.data(), nlat, nlon);
MatrixView<double> vfield(vfield_v.data(), nlat, nlon);
Varray2D<double> urfield(nlevs), vrfield(nlevs);
for (int lid = 0; lid < nlevs; lid++)
......@@ -306,11 +308,15 @@ Mrotuv(void *process)
vrfield[lid].resize(gridsize);
}
Varray<double> uhelp(gridsizex), vhelp(gridsizex);
Varray2D<double> uhelp(nlat, Varray<double>(nlon + 2));
Varray2D<double> vhelp(nlat, Varray<double>(nlon + 2));
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
while (true)
{
const auto nrecs = cdoStreamInqTimestep(streamID1, tsID);
if (nrecs == 0) break;
taxisCopyTimestep(taxisID2, taxisID1);
cdoDefTimestep(streamID2, tsID);
taxisCopyTimestep(taxisID3, taxisID1);
......@@ -318,67 +324,67 @@ Mrotuv(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
int varID, levelID;
cdoInqRecord(streamID1, &varID, &levelID);
if (varID == uid) cdoReadRecord(streamID1, urfield[levelID].data(), &nmiss1);
if (varID == vid) cdoReadRecord(streamID1, vrfield[levelID].data(), &nmiss2);
}
for (levelID = 0; levelID < nlevs; levelID++)
for (int levelID = 0; levelID < nlevs; levelID++)
{
// remove missing values
if (nmiss1 || nmiss2)
{
for (size_t i = 0; i < gridsize; i++)
{
if (DBL_IS_EQUAL(urfield[levelID][i], missval1)) urfield[levelID][i] = 0;
if (DBL_IS_EQUAL(vrfield[levelID][i], missval2)) vrfield[levelID][i] = 0;
if (DBL_IS_EQUAL(urfield[levelID][i], missval1)) urfield[levelID][i] = 0.0;
if (DBL_IS_EQUAL(vrfield[levelID][i], missval2)) vrfield[levelID][i] = 0.0;
}
}
// rotate
rotate_uv(urfield[levelID].data(), vrfield[levelID].data(), nlon, nlat, grid1x.data(), grid1y.data(), ufield.data(),
vfield.data());
rotate_uv(urfield[levelID], vrfield[levelID], nlon, nlat, grid1x, grid1y, ufield_v, vfield_v);
// load to a help field
for (size_t j = 0; j < nlat; j++)
for (size_t i = 0; i < nlon; i++)
{
uhelp[IX2D(j, i + 1, nlon + 2)] = ufield[IX2D(j, i, nlon)];
vhelp[IX2D(j, i + 1, nlon + 2)] = vfield[IX2D(j, i, nlon)];
uhelp[j][i + 1] = ufield[j][i];
vhelp[j][i + 1] = vfield[j][i];
}
// make help field cyclic
for (size_t j = 0; j < nlat; j++)
{
uhelp[IX2D(j, 0, nlon + 2)] = uhelp[IX2D(j, nlon, nlon + 2)];
uhelp[IX2D(j, nlon + 1, nlon + 2)] = uhelp[IX2D(j, 1, nlon + 2)];
vhelp[IX2D(j, 0, nlon + 2)] = vhelp[IX2D(j, nlon, nlon + 2)];
vhelp[IX2D(j, nlon + 1, nlon + 2)] = vhelp[IX2D(j, 1, nlon + 2)];
uhelp[j][0] = uhelp[j][nlon];
uhelp[j][nlon + 1] = uhelp[j][1];
vhelp[j][0] = vhelp[j][nlon];
vhelp[j][nlon + 1] = vhelp[j][1];
}
// interpolate on u/v points
for (size_t j = 0; j < nlat; j++)
for (size_t i = 0; i < nlon; i++)
{
ufield[IX2D(j, i, nlon)] = (uhelp[IX2D(j, i + 1, nlon + 2)] + uhelp[IX2D(j, i + 2, nlon + 2)]) * 0.5;
ufield[j][i] = (uhelp[j][i + 1] + uhelp[j][i + 2]) * 0.5;
}
for (size_t j = 0; j < nlat - 1; j++)
for (size_t i = 0; i < nlon; i++)
{
vfield[IX2D(j, i, nlon)] = (vhelp[IX2D(j, i + 1, nlon + 2)] + vhelp[IX2D(j + 1, i + 1, nlon + 2)]) * 0.5;
vfield[j][i] = (vhelp[j][i + 1] + vhelp[j + 1][i + 1]) * 0.5;
}
for (size_t i = 0; i < nlon; i++)
{
vfield[IX2D(nlat - 1, i, nlon)] = vhelp[IX2D(nlat - 1, i + 1, nlon + 2)];
vfield[nlat - 1][i] = vhelp[nlat - 1][i + 1];
}
cdoDefRecord(streamID2, 0, levelID);
cdoWriteRecord(streamID2, ufield.data(), nmiss1);
cdoWriteRecord(streamID2, ufield_v.data(), nmiss1);
cdoDefRecord(streamID3, 0, levelID);
cdoWriteRecord(streamID3, vfield.data(), nmiss2);
cdoWriteRecord(streamID3, vfield_v.data(), nmiss2);
}
tsID++;
......
......@@ -27,6 +27,7 @@
#include "cdo_options.h"
#include "process_int.h"
#include <mpim_grid.h>
#include "matrix_view.h"
/*
!----------------------------------------------------------------------
......@@ -42,6 +43,7 @@
!
! if this routine is used to rotate data of mpi-om, the
! logical change_sign_v needs to be true.
!
! j. jungclaus: 22.01.04:
! note here for the coupling fields u-i,v_j are on the non-verlapping
! (ie-2=ix) grid, furthermore, the velocity fields were previously
......@@ -50,90 +52,94 @@
*/
void
rotate_uv2(double *u_i, double *v_j, long ix, long iy, double *lon, double *lat, double *u_lon, double *v_lat)
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)
{
/*
in :: u_i(ix,iy),v_j(ix,iy) vector components in i-j-direction
in :: lat(ix,iy),lon(ix,iy) latitudes and longitudes
out :: u_lon(ix,iy),v_lat(ix,iy) vector components in lon-lat direction
in :: u_i[ny][nx], v_j[ny][nx] vector components in i-j-direction
in :: lat[ny][nx], lon[ny][nx] latitudes and longitudes
out :: u_lon[ny][nx], v_lat[ny][nx] vector components in lon-lat direction
*/
double pi = 3.14159265359;
constexpr double pi = 3.14159265359;
MatrixView<double> lon(lon_v.data(), ny, nx);
MatrixView<double> lat(lat_v.data(), ny, nx);
MatrixView<double> u_i(u_i_v.data(), ny, nx);
MatrixView<double> v_j(v_j_v.data(), ny, nx);
MatrixView<double> u_lon(u_lon_v.data(), ny, nx);
MatrixView<double> v_lat(v_lat_v.data(), ny, nx);
// specification whether change in sign is needed for the input arrays
bool change_sign_u = false;
bool change_sign_v = true;
constexpr auto change_sign_u = false;
constexpr auto change_sign_v = true;
// initialization
for (long i = 0; i < ix * iy; i++)
{
v_lat[i] = 0;
u_lon[i] = 0;
}
v_lat_v.assign(nx *ny, 0.0);
u_lon_v.assign(nx *ny, 0.0);
// change sign
if (change_sign_u)
for (long i = 0; i < ix * iy; i++) u_i[i] *= -1;
for (long i = 0; i < nx * ny; i++) u_i_v[i] *= -1;
if (change_sign_v)
for (long i = 0; i < ix * iy; i++) v_j[i] *= -1;
for (long i = 0; i < nx * ny; i++) v_j_v[i] *= -1;
// rotation
for (long j = 0; j < iy; j++)
for (long i = 0; i < ix; i++)
for (long j = 0; j < ny; j++)
for (long i = 0; i < nx; i++)
{
auto ip1 = i + 1;
auto im1 = i - 1;
auto jp1 = j + 1;
auto jm1 = j - 1;
if (ip1 >= ix) ip1 = 0; // the 0-meridian
if (im1 < 0) im1 = ix - 1;
if (jp1 >= iy) jp1 = j; // treatment of the last..
if (ip1 >= nx) ip1 = 0; // the 0-meridian
if (im1 < 0) im1 = nx - 1;
if (jp1 >= ny) jp1 = j; // treatment of the last..
if (jm1 < 0) jm1 = j; // .. and the fist grid-row
// difference in latitudes
auto dlat_i = lat[IX2D(j, ip1, ix)] - lat[IX2D(j, im1, ix)];
auto dlat_j = lat[IX2D(jp1, i, ix)] - lat[IX2D(jm1, i, ix)];
auto dlat_i = lat[j][ip1] - lat[j][im1];
auto dlat_j = lat[jp1][i] - lat[jm1][i];
// difference in longitudes
auto dlon_i = lon[IX2D(j, ip1, ix)] - lon[IX2D(j, im1, ix)];
auto dlon_i = lon[j][ip1] - lon[j][im1];
if (dlon_i > pi) dlon_i -= 2 * pi;
if (dlon_i < (-pi)) dlon_i += 2 * pi;
auto dlon_j = lon[IX2D(jp1, i, ix)] - lon[IX2D(jm1, i, ix)];
auto dlon_j = lon[jp1][i] - lon[jm1][i];
if (dlon_j > pi) dlon_j -= 2 * pi;
if (dlon_j < (-pi)) dlon_j += 2 * pi;
auto lat_factor = std::cos(lat[IX2D(j, i, ix)]);
const auto lat_factor = std::cos(lat[j][i]);
dlon_i = dlon_i * lat_factor;
dlon_j = dlon_j * lat_factor;
// projection by scalar product
u_lon[IX2D(j, i, ix)] = u_i[IX2D(j, i, ix)] * dlon_i + v_j[IX2D(j, i, ix)] * dlat_i;
v_lat[IX2D(j, i, ix)] = u_i[IX2D(j, i, ix)] * dlon_j + v_j[IX2D(j, i, ix)] * dlat_j;
u_lon[j][i] = u_i[j][i] * dlon_i + v_j[j][i] * dlat_i;
v_lat[j][i] = u_i[j][i] * dlon_j + v_j[j][i] * dlat_j;
auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
const auto dist_i = std::sqrt(dlon_i * dlon_i + dlat_i * dlat_i);
const auto dist_j = std::sqrt(dlon_j * dlon_j + dlat_j * dlat_j);
if (std::fabs(dist_i) > 0 && std::fabs(dist_j) > 0)
{
u_lon[IX2D(j, i, ix)] /= dist_i;
v_lat[IX2D(j, i, ix)] /= dist_j;
u_lon[j][i] /= dist_i;
v_lat[j][i] /= dist_j;
}
else
{
u_lon[IX2D(j, i, ix)] = 0.0;
v_lat[IX2D(j, i, ix)] = 0.0;
u_lon[j][i] = 0.0;
v_lat[j][i] = 0.0;
}
if (Options::cdoVerbose)
{
// velocity vector lengths
auto absold = std::sqrt(u_i[IX2D(j, i, ix)] * u_i[IX2D(j, i, ix)] + v_j[IX2D(j, i, ix)] * v_j[IX2D(j, i, ix)]);
auto absnew = std::sqrt(u_lon[IX2D(j, i, ix)] * u_lon[IX2D(j, i, ix)] + v_lat[IX2D(j, i, ix)] * v_lat[IX2D(j, i, ix)]);
const auto absold = std::sqrt(u_i[j][i] * u_i[j][i] + v_j[j][i] * v_j[j][i]);
const auto absnew = std::sqrt(u_lon[j][i] * u_lon[j][i] + v_lat[j][i] * v_lat[j][i]);
if (i % 20 ==