Commit f19494af authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Vertint*: changed type of vert_index from int to long.

parent a246702b
Pipeline #4150 passed with stages
in 14 minutes and 38 seconds
......@@ -3,6 +3,14 @@
* Using CDI library version 1.9.9
* Version 1.9.9 release
2020-09-10 Uwe Schulzweida
* Vertint*: changed type of vert_index from int to long
2020-09-09 Uwe Schulzweida
* Set MAX_PROCESS to 512
2020-08-26 Uwe Schulzweida
* setgridtype,regular: bug fix for regional reduced Gaussian grids
......
......@@ -257,7 +257,7 @@ Vertintap(void *process)
if (nlev != nhlev) cdoAbort("Internal error, wrong number of height level!");
}
std::vector<int> vert_index;
std::vector<long> vert_index;
Field ps_prog;
Field3D full_press, half_press;
if (zaxisIDh != -1 && gridsize > 0)
......@@ -425,6 +425,7 @@ Vertintap(void *process)
}
vertical_interp_X(nlevels, full_press, half_press, vardata1[varID], vardata2[varID], vert_index, plev, gridsize);
if (!extrapolate) varrayCopy(nplev, pnmiss, varnmiss[varID]);
}
......
......@@ -233,7 +233,6 @@ is_height_axis(int zaxisID)
void *
Vertintgh(void *process)
{
#define NEWIMPL 1
int nrecs;
int varID, levelID;
int f_heightID = -1, h_heightID = -1;;
......@@ -357,38 +356,14 @@ Vertintgh(void *process)
varListInit(varList2, vlistID2);
varListSetMemtype(varList2, memtype);
#ifdef NEWIMPL
Field height_bottom;
std::vector<size_t> pnmiss;
if (!extrapolate) pnmiss.resize(nheightlevs);
std::vector<int> vert_index_full, vert_index_half;
std::vector<long> vert_index_full, vert_index_half;
if (-1 != f_heightID) vert_index_full.resize(gridsize * nheightlevs);
if (-1 != h_heightID) vert_index_half.resize(gridsize * nheightlevs);
#else
Varray<int> lev_idx1, lev_idx2;
Varray<double> lev_wgt1, lev_wgt2;
if (-1 != f_heightID)
{
const auto wsize = gridsize * nfulllevs;
lev_idx1.resize(wsize);
lev_idx2.resize(wsize);
lev_wgt1.resize(wsize);
lev_wgt2.resize(wsize);
}
Varray<int> xlev_idx1, xlev_idx2;
Varray<double> xlev_wgt1, xlev_wgt2;
if (-1 != h_heightID)
{
const auto wsize = gridsize * nhalflevs;
xlev_idx1.resize(wsize);
xlev_idx2.resize(wsize);
xlev_wgt1.resize(wsize);
xlev_wgt2.resize(wsize);
}
#endif
std::vector<bool> vars(nvars), varinterp(nvars);
std::vector<std::vector<size_t>> varnmiss(nvars);
......@@ -461,7 +436,6 @@ Vertintgh(void *process)
for (varID = 0; varID < nvars; varID++)
if (varinterp[varID]) vars[varID] = true;
#ifdef NEWIMPL
bool lreverse = true;
if (-1 != f_heightID && (tsID == 0 || varList1[f_heightID].timetype != TIME_CONSTANT))
{
......@@ -473,6 +447,7 @@ Vertintgh(void *process)
genindmiss(vert_index_full, heightlevs, gridsize, height_bottom, pnmiss, lreverse);
}
}
if (-1 != h_heightID && (tsID == 0 || varList1[h_heightID].timetype != TIME_CONSTANT))
{
genind(vert_index_half, heightlevs, vardata1[h_heightID], gridsize, lreverse);
......@@ -483,13 +458,6 @@ Vertintgh(void *process)
genindmiss(vert_index_half, heightlevs, gridsize, height_bottom, pnmiss, lreverse);
}
}
#else
if (-1 != f_heightID && (tsID == 0 || varList1[f_heightID].timetype != TIME_CONSTANT))
vert_gen_weights3d1d(extrapolate, gridsize, nfulllevs, vardata1[f_heightID], nheightlevs, heightlevs, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
if (-1 != h_heightID && (tsID == 0 || varList1[h_heightID].timetype != TIME_CONSTANT))
vert_gen_weights3d1d(extrapolate, gridsize, nhalflevs, vardata1[h_heightID], nheightlevs, heightlevs, xlev_idx1, xlev_idx2, xlev_wgt1, xlev_wgt2);
#endif
for (varID = 0; varID < nvars; varID++)
{
......@@ -508,30 +476,12 @@ Vertintgh(void *process)
if (varnmiss[varID][levelID]) cdoAbort("Missing values unsupported for this operator!");
}
#ifdef NEWIMPL
if (nlevels == nfulllevs)
vertical_interp_X(nlevels, vardata1[f_heightID], vardata1[varID], vardata2[varID], vert_index_full, heightlevs, gridsize);
else
vertical_interp_X(nlevels, vardata1[h_heightID], vardata1[varID], vardata2[varID], vert_index_half, heightlevs, gridsize);
if (!extrapolate) varrayCopy(nheightlevs, pnmiss, varnmiss[varID]);
#else
const auto missval = varList1[varID].missval;
if (nlevels == nfulllevs)
vert_interp_lev3d(gridsize, missval, vardata1[varID], vardata2[varID], nheightlevs, lev_idx1, lev_idx2, lev_wgt1, lev_wgt2);
else
vert_interp_lev3d(gridsize, missval, vardata1[varID], vardata2[varID], nheightlevs, xlev_idx1, xlev_idx2, xlev_wgt1, xlev_wgt2);
if (!extrapolate)
for (levelID = 0; levelID < nheightlevs; levelID++)
{
const auto offset = gridsize * levelID;
if (memtype == MemType::Float)
varnmiss[varID][levelID] = arrayNumMV(gridsize, &vardata2[varID].vec_f[offset], (float)missval);
else
varnmiss[varID][levelID] = arrayNumMV(gridsize, &vardata2[varID].vec_d[offset], missval);
}
#endif
}
for (levelID = 0; levelID < varList2[varID].nlevels; levelID++)
......
......@@ -227,7 +227,7 @@ Vertintml(void *process)
if (nlev != nhlev) cdoAbort("Internal error, wrong number of hybrid level!");
}
std::vector<int> vert_index;
std::vector<long> vert_index;
Field3D full_press, half_press;
if (zaxisIDh != -1 && gridsize > 0)
......
......@@ -83,7 +83,7 @@ struct Control
int nvct;
double *vct;
int *vert_index;
long *vert_index;
size_t *pnmiss;
double *Orography;
double *p_of_height;
......
......@@ -1808,7 +1808,7 @@ after_processML(struct Control *globs, struct Variable *vars)
if (globs->Type >= 30)
{
if (globs->vert_index == nullptr) globs->vert_index = (int *) Malloc(globs->NumLevelRequest * globs->DimGP * sizeof(int));
if (globs->vert_index == nullptr) globs->vert_index = (long *) Malloc(globs->NumLevelRequest * globs->DimGP * sizeof(long));
if (globs->unitsel)
{
......
......@@ -18,7 +18,7 @@
#include "field_vinterp.h"
void
genind(std::vector<int> &vert_index, Varray<double> &plev, Field3D &full_level, size_t gridsize, bool lreverse)
genind(std::vector<long> &vert_index, Varray<double> &plev, Field3D &full_level, size_t gridsize, bool lreverse)
{
const auto nplev = plev.size();
const auto nhlevf = full_level.nlevels;
......@@ -29,7 +29,7 @@ genind(std::vector<int> &vert_index, Varray<double> &plev, Field3D &full_level,
}
void
genindmiss(std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize, Field &level0, std::vector<size_t> &pnmiss, bool lreverse)
genindmiss(std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize, Field &level0, std::vector<size_t> &pnmiss, bool lreverse)
{
const auto nplev = plev.size();
if (level0.memType == MemType::Float)
......@@ -39,7 +39,7 @@ genindmiss(std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize,
}
void
vertical_interp_T(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field &sgeopot, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize)
vertical_interp_T(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field &sgeopot, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize)
{
const auto nplev = plev.size();
const auto missval = field1.missval;
......@@ -56,7 +56,7 @@ vertical_interp_T(size_t nlevels, Field3D &full_level, Field3D &half_level, Fiel
}
void
vertical_interp_Z(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field3D &temp, Field &sgeopot, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize)
vertical_interp_Z(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field3D &temp, Field &sgeopot, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize)
{
const auto nplev = plev.size();
const auto missval = field1.missval;
......@@ -73,7 +73,7 @@ vertical_interp_Z(size_t nlevels, Field3D &full_level, Field3D &half_level, Fiel
}
void
vertical_interp_X(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize)
vertical_interp_X(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize)
{
const auto nplev = plev.size();
const auto missval = field1.missval;
......@@ -92,7 +92,7 @@ vertical_interp_X(size_t nlevels, Field3D &full_level, Field3D &half_level, Fiel
}
void
vertical_interp_X(size_t nlevels, Field3D &levels, Field3D &field1, Field3D &field2, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize)
vertical_interp_X(size_t nlevels, Field3D &levels, Field3D &field1, Field3D &field2, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize)
{
const auto nplev = plev.size();
const auto missval = field1.missval;
......
......@@ -21,16 +21,16 @@
#include "field.h"
#include "vertical_interp.h"
void genind(std::vector<int> &vert_index, Varray<double> &plev, Field3D &full_level, size_t gridsize, bool lreverse = false);
void genind(std::vector<long> &vert_index, Varray<double> &plev, Field3D &full_level, size_t gridsize, bool lreverse = false);
void genindmiss(std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize, Field &level0, std::vector<size_t> &pnmiss, bool lreverse = false);
void genindmiss(std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize, Field &level0, std::vector<size_t> &pnmiss, bool lreverse = false);
void vertical_interp_T(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field &sgeopot, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_T(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field &sgeopot, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_Z(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field3D &temp, Field &sgeopot, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_Z(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, Field3D &temp, Field &sgeopot, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_X(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_X(size_t nlevels, Field3D &full_level, Field3D &half_level, Field3D &field1, Field3D &field2, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_X(size_t nlevels, Field3D &levels, Field3D &field1, Field3D &field2, std::vector<int> &vert_index, Varray<double> &plev, size_t gridsize);
void vertical_interp_X(size_t nlevels, Field3D &levels, Field3D &field1, Field3D &field2, std::vector<long> &vert_index, Varray<double> &plev, size_t gridsize);
#endif
......@@ -204,7 +204,7 @@ extrapolate_Z(double pres, double halfp, double fullp, double geop, double temp)
template <typename T>
void
vertical_interp_T(const T *restrict geop, const T *restrict gt, T *pt, const T *restrict fullp,
const T *restrict halfp, const int *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
const T *restrict halfp, const long *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
double missval)
{
#ifdef _OPENMP
......@@ -214,7 +214,7 @@ vertical_interp_T(const T *restrict geop, const T *restrict gt, T *pt, const T *
{
long nl, nh;
const auto pres = plev[lp];
const int *restrict nxl = nx + lp * ngp;
const long *restrict nxl = nx + lp * ngp;
auto ptl = pt + lp * ngp;
for (long i = 0; i < ngp; i++)
......@@ -242,16 +242,16 @@ vertical_interp_T(const T *restrict geop, const T *restrict gt, T *pt, const T *
}
template void vertical_interp_T(const float *restrict geop, const float *restrict gt, float *pt, const float *restrict fullp,
const float *restrict halfp, const int *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
const float *restrict halfp, const long *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
double missval);
template void vertical_interp_T(const double *restrict geop, const double *restrict gt, double *pt, const double *restrict fullp,
const double *restrict halfp, const int *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
const double *restrict halfp, const long *nx, const double *restrict plev, long nplev, long ngp, long nhlev,
double missval);
template <typename T>
void
vertical_interp_Z(const T *restrict geop, const T *restrict gz, T *pz, const T *restrict fullp,
const T *restrict halfp, const int *nx, const T *restrict gt, const double *restrict plev, long nplev,
const T *restrict halfp, const long *nx, const T *restrict gt, const double *restrict plev, long nplev,
long ngp, long nhlev, double missval)
{
assert(geop != nullptr);
......@@ -267,7 +267,7 @@ vertical_interp_Z(const T *restrict geop, const T *restrict gz, T *pz, const T *
{
long nl, nh;
const auto pres = plev[lp];
const int *restrict nxl = nx + lp * ngp;
const long *restrict nxl = nx + lp * ngp;
auto pzl = pz + lp * ngp;
for (long i = 0; i < ngp; i++)
......@@ -297,10 +297,10 @@ vertical_interp_Z(const T *restrict geop, const T *restrict gz, T *pz, const T *
}
template void vertical_interp_Z(const float *restrict geop, const float *restrict gz, float *pz, const float *restrict fullp,
const float *restrict halfp, const int *nx, const float *restrict gt, const double *restrict plev, long nplev,
const float *restrict halfp, const long *nx, const float *restrict gt, const double *restrict plev, long nplev,
long ngp, long nhlev, double missval);
template void vertical_interp_Z(const double *restrict geop, const double *restrict gz, double *pz, const double *restrict fullp,
const double *restrict halfp, const int *nx, const double *restrict gt, const double *restrict plev, long nplev,
const double *restrict halfp, const long *nx, const double *restrict gt, const double *restrict plev, long nplev,
long ngp, long nhlev, double missval);
template <typename T>
......@@ -316,7 +316,7 @@ vertical_interp_X_kernel(const T *restrict gt, const T *restrict hyb_press, long
template <typename T>
void
vertical_interp_X(const T *restrict gt, T *pt, const T *hyb_press, const int *nx, const double *restrict plev,
vertical_interp_X(const T *restrict gt, T *pt, const T *hyb_press, const long *nx, const double *restrict plev,
long nplev, long ngp, long nhlev, double missval)
{
if (nplev > 3)
......@@ -327,7 +327,7 @@ vertical_interp_X(const T *restrict gt, T *pt, const T *hyb_press, const int *nx
for (long lp = 0; lp < nplev; lp++)
{
auto pres = plev[lp];
const int *restrict nxl = nx + lp * ngp;
const long *restrict nxl = nx + lp * ngp;
auto ptl = pt + lp * ngp;
for (long i = 0; i < ngp; i++)
{
......@@ -340,7 +340,7 @@ vertical_interp_X(const T *restrict gt, T *pt, const T *hyb_press, const int *nx
for (long lp = 0; lp < nplev; lp++)
{
auto pres = plev[lp];
const int *restrict nxl = nx + lp * ngp;
const long *restrict nxl = nx + lp * ngp;
auto *restrict ptl = pt + lp * ngp;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(gt, ptl, hyb_press, nxl, pres, ngp, nhlev, missval)
......@@ -353,17 +353,17 @@ vertical_interp_X(const T *restrict gt, T *pt, const T *hyb_press, const int *nx
}
}
template void vertical_interp_X(const float *restrict gt, float *pt, const float *hyb_press, const int *nx,
template void vertical_interp_X(const float *restrict gt, float *pt, const float *hyb_press, const long *nx,
const double *restrict plev, long nplev, long ngp, long nhlev, double missval);
template void vertical_interp_X(const double *restrict gt, double *pt, const double *hyb_press, const int *nx,
template void vertical_interp_X(const double *restrict gt, double *pt, const double *hyb_press, const long *nx,
const double *restrict plev, long nplev, long ngp, long nhlev, double missval);
template <typename T>
void
genind(int *nx, const double *restrict plev, const T *restrict fullp, long ngp, long nplev, long nhlev, bool lreverse)
genind(long *nx, const double *restrict plev, const T *restrict fullp, long ngp, long nplev, long nhlev, bool lreverse)
{
varrayFill(ngp * nplev, nx, 0);
varrayFill(ngp * nplev, nx, 0L);
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(nx, plev, fullp, ngp, nplev, nhlev, lreverse)
......@@ -371,7 +371,7 @@ genind(int *nx, const double *restrict plev, const T *restrict fullp, long ngp,
for (long lp = 0; lp < nplev; lp++)
{
const T pres = plev[lp];
int *restrict nxl = nx + lp * ngp;
long *restrict nxl = nx + lp * ngp;
for (long lh = 0; lh < nhlev; lh++)
{
const auto *restrict fullpx = fullp + lh * ngp;
......@@ -394,12 +394,12 @@ genind(int *nx, const double *restrict plev, const T *restrict fullp, long ngp,
}
// Explicit instantiation
template void genind(int *nx, const double *plev, const float *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
template void genind(int *nx, const double *plev, const double *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
template void genind(long *nx, const double *plev, const float *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
template void genind(long *nx, const double *plev, const double *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
template <typename T>
void
genindmiss(int *nx, const double *restrict plev, int ngp, int nplev, const T *restrict ps_prog, size_t *restrict pnmiss, bool lreverse)
genindmiss(long *nx, const double *restrict plev, long ngp, long nplev, const T *restrict ps_prog, size_t *restrict pnmiss, bool lreverse)
{
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(nx, plev, ngp, nplev, ps_prog, pnmiss, lreverse)
......@@ -408,7 +408,7 @@ genindmiss(int *nx, const double *restrict plev, int ngp, int nplev, const T *re
{
pnmiss[lp] = 0;
const T pres = plev[lp];
int *restrict nxl = nx + lp * ngp;
long *restrict nxl = nx + lp * ngp;
if (lreverse)
{
......@@ -437,5 +437,5 @@ genindmiss(int *nx, const double *restrict plev, int ngp, int nplev, const T *re
}
// Explicit instantiation
template void genindmiss(int *nx, const double *plev, int ngp, int nplev, const float *ps_prog, size_t *pnmiss, bool lreverse = false);
template void genindmiss(int *nx, const double *plev, int ngp, int nplev, const double *ps_prog, size_t *pnmiss, bool lreverse = false);
template void genindmiss(long *nx, const double *plev, long ngp, long nplev, const float *ps_prog, size_t *pnmiss, bool lreverse = false);
template void genindmiss(long *nx, const double *plev, long ngp, long nplev, const double *ps_prog, size_t *pnmiss, bool lreverse = false);
......@@ -27,21 +27,21 @@ void presh(T *fullp, T *halfp, const double *vct, const T *ps, long nhlev, long
void extrapolate_P(double *slp, const double *halfp, const double *fullp, const double *geop, const double *temp, long ngp);
template <typename T>
void vertical_interp_T(const T *geop, const T *gt, T *pt, const T *fullp, const T *halfp, const int *nx,
void vertical_interp_T(const T *geop, const T *gt, T *pt, const T *fullp, const T *halfp, const long *nx,
const double *plev, long nplev, long ngp, long nhlev, double missval);
template <typename T>
void vertical_interp_Z(const T *geop, const T *gz, T *pz, const T *fullp, const T *halfp, const int *nx,
void vertical_interp_Z(const T *geop, const T *gz, T *pz, const T *fullp, const T *halfp, const long *nx,
const T *gt, const double *plev, long nplev, long ngp, long nhlev, double missval);
template <typename T>
void vertical_interp_X(const T *gt, T *pt, const T *hyb_press, const int *nx, const double *plev,
void vertical_interp_X(const T *gt, T *pt, const T *hyb_press, const long *nx, const double *plev,
long nplev, long ngp, long nhlev, double missval);
template <typename T>
void genind(int *nx, const double *plev, const T *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
void genind(long *nx, const double *plev, const T *fullp, long ngp, long nplev, long nhlev, bool lreverse = false);
template <typename T>
void genindmiss(int *nx, const double *plev, int ngp, int nplev, const T *ps_prog, size_t *pnmiss, bool lreverse = false);
void genindmiss(long *nx, const double *plev, long ngp, long nplev, const T *ps_prog, size_t *pnmiss, bool lreverse = false);
#endif /* VERTICAL_INTERP_H */
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