Commit 86a77401 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

vert_gen_weights3d: OpenMP version.

parent edb297e7
Pipeline #4340 failed with stages
in 13 minutes and 22 seconds
......@@ -24,6 +24,7 @@
#include <cdi.h>
#include "cimdOmp.h"
#include "cdo_options.h"
#include "process_int.h"
#include "cdo_vlist.h"
......@@ -48,32 +49,38 @@ static void
vert_gen_weights3d(bool expol, size_t gridsize, int nlev1, Varray<float> &xlev1, int nlev2, Varray<double> &xlev2,
Varray<int> &xlev_idx1, Varray<int> &xlev_idx2, Varray<float> &xlev_wgt1, Varray<float> &xlev_wgt2)
{
Varray<double> lev1p2(nlev1 + 2);
Varray<double> lev2(nlev2);
Varray<float> lev_wgt1(nlev2), lev_wgt2(nlev2);
Varray<int> lev_idx1(nlev2), lev_idx2(nlev2);
const auto nthreads = Threading::ompNumThreads;
Varray2D<double> lev1p2(nthreads, Varray<double>(nlev1 + 2));
Varray2D<double> lev2(nthreads, Varray<double>(nlev2));
Varray2D<float> lev_wgt1(nthreads, Varray<float>(nlev2)), lev_wgt2(nthreads, Varray<float>(nlev2));
Varray2D<int> lev_idx1(nthreads, Varray<int>(nlev2)), lev_idx2(nthreads, Varray<int>(nlev2));
// Check monotony of vertical levels
for (int k = 0; k < nlev1; ++k) lev1p2[k] = xlev1[k * gridsize];
const auto lup = levelDirUp(nlev1, lev1p2.data());
const auto ldown = levelDirDown(nlev1, lev1p2.data());
for (int k = 0; k < nlev1; ++k) lev1p2[0][k] = xlev1[k * gridsize];
const auto lup = levelDirUp(nlev1, lev1p2[0].data());
const auto ldown = levelDirDown(nlev1, lev1p2[0].data());
if (!lup && !ldown) cdoAbort("Non monotonic zaxis!");
const double level_0 = lup ? -1.e33 : 1.e33;
const double level_N = lup ? 1.e33 : -1.e33;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(expol, gridsize, nlev1, nlev2, lev2, xlev2, level_0, level_N, lev1p2, xlev1, xlev_idx1, xlev_idx2, xlev_wgt1, xlev_wgt2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2)
#endif
for (size_t i = 0; i < gridsize; i++)
{
lev1p2[0] = level_0;
lev1p2[nlev1 + 1] = level_N;
for (int k = 0; k < nlev1; ++k) lev1p2[k + 1] = xlev1[k * gridsize + i];
for (int k = 0; k < nlev2; ++k) lev2[k] = xlev2[k * gridsize + i];
const auto ompthID = cdo_omp_get_thread_num();
vert_gen_weights(expol, nlev1 + 2, lev1p2, nlev2, lev2, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
lev1p2[ompthID][0] = level_0;
lev1p2[ompthID][nlev1 + 1] = level_N;
for (int k = 0; k < nlev1; ++k) lev1p2[ompthID][k + 1] = xlev1[k * gridsize + i];
for (int k = 0; k < nlev2; ++k) lev2[ompthID][k] = xlev2[k * gridsize + i];
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[k];
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[k];
for (int k = 0; k < nlev2; ++k) xlev_wgt1[k * gridsize + i] = lev_wgt1[k];
for (int k = 0; k < nlev2; ++k) xlev_wgt2[k * gridsize + i] = lev_wgt2[k];
vert_gen_weights(expol, nlev1 + 2, lev1p2[ompthID], nlev2, lev2[ompthID], lev_idx1[ompthID], lev_idx2[ompthID], lev_wgt1[ompthID], lev_wgt2[ompthID]);
for (int k = 0; k < nlev2; ++k) xlev_idx1[k * gridsize + i] = lev_idx1[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_idx2[k * gridsize + i] = lev_idx2[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_wgt1[k * gridsize + i] = lev_wgt1[ompthID][k];
for (int k = 0; k < nlev2; ++k) xlev_wgt2[k * gridsize + i] = lev_wgt2[ompthID][k];
}
}
......@@ -100,7 +107,7 @@ Intlevel3d(void *process)
const bool expol = (operatorID == INTLEVELX3D);
operatorInputArg("icoordinate");
operatorInputArg("tgtcoordinate");
const auto streamID1 = cdoOpenRead(0); // input data
const auto streamID3 = cdoOpenWrite(2); // output stream
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment