Commit 30e904dc authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added class FC_Transformation.

parent 9adf398d
......@@ -249,10 +249,10 @@ case "${HOSTNAME}" in
# stretch workstation x64
breeze*)
CDOLIBS="--with-netcdf=/sw/stretch-x64/netcdf/netcdf_c-4.6.1 \
--with-fftw3 \
--with-udunits2=/sw/stretch-x64/sys/udunits-2.2.24-gccsys"
if test "$COMP" = intel ; then
${CONFPATH}configure \
--with-fftw3 \
$CDOLIBS \
CC=icc CXX=icpc CFLAGS="-g -Wall -O2 -qopt-report=5 -march=native"
elif test "$COMP" = pgi ; then
......@@ -284,7 +284,6 @@ case "${HOSTNAME}" in
CC=gcc CFLAGS="-g -O2 -fstack-protector-strong"
else
${CONFPATH}configure \
--with-fftw3 \
$CDOLIBS \
CXX=g++ CXXFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native" \
CC=gcc CFLAGS="-g -Wall -W -Wfloat-equal -pedantic -O3 -march=native"
......@@ -297,6 +296,7 @@ case "${HOSTNAME}" in
--with-hdf5=/sw/rhel6-x64/hdf5/hdf5-1.8.14-threadsafe-gcc48 \
--with-szlib=/sw/rhel6-x64/sys/libaec-0.3.4-gcc48 \
--with-udunits2=/sw/rhel6-x64/util/udunits-2.2.17-gcc48 \
--with-fftw3 \
--with-curl \
--with-ossp-uuid \
--with-proj=/sw/rhel6-x64/graphics/proj5-5.2.0-gcc64 \
......@@ -305,7 +305,6 @@ case "${HOSTNAME}" in
if test "$COMP" = intel ; then
${CONFPATH}configure --prefix=$HOME/local \
--with-fftw3 \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,/sw/rhel6-x64/eccodes/eccodes-2.6.0-gcc64/lib" \
F77=ifort FFLAGS="-g -O2" \
......@@ -319,7 +318,6 @@ case "${HOSTNAME}" in
CC=pgcc CFLAGS="-g -fast"
elif test "$COMP" = gnu_pic ; then
${CONFPATH}configure --prefix=$HOME/local \
--with-fftw3 \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,/sw/rhel6-x64/eccodes/eccodes-2.6.0-gcc64/lib" \
F77=gfortran FFLAGS="-g -O2" \
......@@ -327,7 +325,6 @@ case "${HOSTNAME}" in
CC=gcc CFLAGS='-g -Wall -O3 -march=native -mavx2 -fPIC'
else
${CONFPATH}configure --prefix=$HOME/local \
--with-fftw3 \
$CDOLIBS \
LDFLAGS="-Wl,-rpath,/sw/rhel6-x64/eccodes/eccodes-2.6.0-gcc64/lib" \
F77=gfortran FFLAGS="-g -O2" \
......
......@@ -180,8 +180,15 @@ FC(void *process)
}
}
SPTRANS sptrans;
if (nlon > 0) sptrans.init(nlon, nlat, ntr, 0);
FC_Transformation fcTrans;
SP_Transformation spTrans;
if (nlon > 0)
{
if (operatorID == GP2FC || operatorID == FC2GP)
fcTrans.init(nlon, nlat, ntr);
else
spTrans.init(nlon, nlat, ntr, 0);
}
// printf("nfc %d, ntr %d, nlat %zu, nlon %zu\n", nfc, ntr, nlat, nlon);
......@@ -221,13 +228,13 @@ FC(void *process)
gridID1 = vlistInqVarGrid(vlistID1, varID);
if (operatorID == FC2SP)
four2spec(sptrans, gridID1, array1.data(), gridID2, array2.data());
four2spec(spTrans, gridID1, array1.data(), gridID2, array2.data());
else if (operatorID == SP2FC)
spec2four(sptrans, gridID1, array1.data(), gridID2, array2.data());
spec2four(spTrans, gridID1, array1.data(), gridID2, array2.data());
else if (operatorID == FC2GP)
four2grid(sptrans, gridID1, array1.data(), gridID2, array2.data());
four2grid(fcTrans, gridID1, array1.data(), gridID2, array2.data());
else if (operatorID == GP2FC)
grid2four(sptrans, gridID1, array1.data(), gridID2, array2.data());
grid2four(fcTrans, gridID1, array1.data(), gridID2, array2.data());
cdoDefRecord(streamID2, varID, levelID);
cdoWriteRecord(streamID2, array2.data(), nmiss);
......
......@@ -45,7 +45,7 @@ Spectral(void *process)
int gridID1 = -1, gridID2 = -1;
size_t nmiss;
std::vector<int> waves;
SPTRANS sptrans;
SP_Transformation spTrans;
ListArray<int> listArrayInt;
cdoInitialize(process);
......@@ -123,7 +123,7 @@ Spectral(void *process)
gridID2 = gridIDsp;
ntr = gridInqTrunc(gridID2);
sptrans.init(nlon, nlat, ntr, 0);
spTrans.init(nlon, nlat, ntr, 0);
}
}
else if (lsp2gp)
......@@ -153,7 +153,7 @@ Spectral(void *process)
const long ntr = gridInqTrunc(gridID1);
const long nlon = gridInqXsize(gridID2);
const long nlat = gridInqYsize(gridID2);
sptrans.init(nlon, nlat, ntr, 0);
spTrans.init(nlon, nlat, ntr, 0);
}
}
else if (operatorID == SP2SP)
......@@ -230,9 +230,9 @@ Spectral(void *process)
gridID1 = vlistInqVarGrid(vlistID1, varID);
if (lgp2sp)
grid2spec(sptrans, gridID1, array1.data(), gridID2, array2.data());
grid2spec(spTrans, gridID1, array1.data(), gridID2, array2.data());
else if (lsp2gp)
spec2grid(sptrans, gridID1, array1.data(), gridID2, array2.data());
spec2grid(spTrans, gridID1, array1.data(), gridID2, array2.data());
else if (operatorID == SP2SP)
spec2spec(gridID1, array1.data(), gridID2, array2.data());
else if (operatorID == SPCUT)
......
......@@ -45,8 +45,8 @@ Wind(void *process)
size_t nmiss;
long ntr = -1;
int varID1 = -1, varID2 = -1;
SPTRANS sptrans;
DVTRANS dvtrans;
SP_Transformation spTrans;
DV_Transformation dvTrans;
char varname[CDI_MAX_NAME];
cdoInitialize(process);
......@@ -208,7 +208,7 @@ Wind(void *process)
const long nlat = gridInqYsize(gridID1);
ntr = gridInqTrunc(gridID2);
sptrans.init(nlon, nlat, ntr, 1);
spTrans.init(nlon, nlat, ntr, 1);
}
}
else if (ldv2uv)
......@@ -256,8 +256,8 @@ Wind(void *process)
const long nlat = gridInqYsize(gridID2);
ntr = gridInqTrunc(gridID1);
sptrans.init(nlon, nlat, ntr, 0);
dvtrans.init(ntr);
spTrans.init(nlon, nlat, ntr, 0);
dvTrans.init(ntr);
}
}
else if (operatorID == DV2PS)
......@@ -349,9 +349,9 @@ Wind(void *process)
if (varID1 != -1 && varID2 != -1)
{
if (luv2dv)
trans_uv2dv(sptrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
trans_uv2dv(spTrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
else if (ldv2uv)
trans_dv2uv(sptrans, dvtrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
trans_dv2uv(spTrans, dvTrans, nlev, gridID1, ivar1.data(), ivar2.data(), gridID2, ovar1.data(), ovar2.data());
else if (operatorID == DV2PS)
{
dv2ps(ivar1.data(), ovar1.data(), nlev, ntr);
......
......@@ -129,7 +129,7 @@ dv2uv_kernel(const double *d, const double *o, double *restrict u, double *restr
}
void
dv2uv(const double *d, const double *o, double *u, double *v, double *f, double *g, long nt, long nsp, long nlev)
dv2uv(const double *d, const double *o, double *u, double *v, const double *f, const double *g, long nt, long nsp, long nlev)
{
// d(nsp,nlev), o(nsp,nlev) ! divergence, vorticity
// u(nsp,nlev), v(nsp,nlev) ! zonal wind, meridional wind
......@@ -176,42 +176,35 @@ scaluv(double *fu, const double *rclat, long nlat, long lot)
}
void
uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, long klev, long nlat, long nt)
uv2dv(const double *fu, const double *fv, double *sd, double *sv, const double *pol2, const double *pol3, long klev, long nlat, long nt)
{
long lev, jmm, jfc, lat;
double dir, dii, vor, voi;
double *ufr, *ufi, *vfr, *vfi;
double *ful, *fvl, *sdl, *svl;
double *po2, *po3;
long nsp2 = (nt + 1) * (nt + 2);
long nfc = (nt + 1) * 2;
#if defined(_OPENMP)
#pragma omp parallel for default(shared) private(jmm, jfc, lat, po2, po3, ful, fvl, sdl, svl, ufr, ufi, vfr, vfi, dir, dii, vor, \
voi)
#pragma omp parallel for default(shared)
#endif
for (lev = 0; lev < klev; lev++)
for (long lev = 0; lev < klev; lev++)
{
po2 = pol2;
po3 = pol3;
ful = fu + lev * nfc * nlat;
fvl = fv + lev * nfc * nlat;
sdl = sd + lev * nsp2;
svl = sv + lev * nsp2;
for (jmm = 0; jmm <= nt; jmm++)
auto po2 = pol2;
auto po3 = pol3;
auto ful = fu + lev * nfc * nlat;
auto fvl = fv + lev * nfc * nlat;
auto sdl = sd + lev * nsp2;
auto svl = sv + lev * nsp2;
for (long jmm = 0; jmm <= nt; jmm++)
{
for (jfc = jmm; jfc <= nt; jfc++)
for (long jfc = jmm; jfc <= nt; jfc++)
{
ufr = ful;
ufi = ful + nlat;
vfr = fvl;
vfi = fvl + nlat;
dir = 0.0;
dii = 0.0;
vor = 0.0;
voi = 0.0;
for (lat = 0; lat < nlat; lat++)
auto ufr = ful;
auto ufi = ful + nlat;
auto vfr = fvl;
auto vfi = fvl + nlat;
double dir = 0.0;
double dii = 0.0;
double vor = 0.0;
double voi = 0.0;
for (long lat = 0; lat < nlat; lat++)
{
dir += vfr[lat] * po2[lat] - ufi[lat] * po3[lat];
dii += vfi[lat] * po2[lat] + ufr[lat] * po3[lat];
......@@ -234,15 +227,12 @@ uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3
void
geninx(long ntr, double *f, double *g)
{
long m2, n2;
long m, n;
for (m = 0; m <= ntr; m++)
for (long m = 0; m <= ntr; m++)
{
m2 = m * m;
for (n = m; n <= ntr; n++)
long m2 = m * m;
for (long n = m; n <= ntr; n++)
{
n2 = n * n;
long n2 = n * n;
if (n)
{
*g++ = -PlanetRadius / n * std::sqrt((double) (n2 - m2) / (double) (4 * n2 - 1));
......
......@@ -1910,7 +1910,7 @@ qpassc(const double *a, const double *b, double *c, double *d, const double *tri
/* Fast Fourier Transform */
/* ====================== */
void
fc2gp(double *restrict trig, const long *restrict ifax, double *restrict fc, double *restrict gp, long nlat, long nlon, long nlev,
fc2gp(const double *restrict trig, const long *restrict ifax, double *restrict fc, double *restrict gp, long nlat, long nlon, long nlev,
long nfc)
{
/* fc2gp performs fourier to gridpoint transforms using */
......@@ -2057,7 +2057,7 @@ fc2gp(double *restrict trig, const long *restrict ifax, double *restrict fc, dou
}
void
gp2fc(double *trig, const long *ifax, const double *restrict gp, double *restrict fc, long nlat, long nlon, long nlev, long nfc)
gp2fc(const double *trig, const long *ifax, const double *restrict gp, double *restrict fc, long nlat, long nlon, long nlev, long nfc)
{
long lot, fou, ia, ifac, jump, k, la;
long lat, lev, lon, nfax, rix, wix;
......
......@@ -359,8 +359,8 @@ after_legini(long ntr, long nlat, double *restrict poli, double *restrict pold,
{
#ifdef _OPENMP
const auto ompthID = cdo_omp_get_thread_num();
double *pnm = pnm2[ompthID].data();
double *work = work2[ompthID].data();
auto pnm = pnm2[ompthID].data();
auto work = work2[ompthID].data();
#endif
const double zgwt = gwt[jgl];
......@@ -399,9 +399,9 @@ sp2fctest(const double *sa, double *fa, const double *poli, long nlev, long nlat
for (long lev = 0; lev < nlev; lev++)
{
const double *restrict pol = poli;
const double *restrict sal = sa + lev * nsp2;
double *fal = fa + lev * nfc * nlat;
auto pol = poli;
auto sal = sa + lev * nsp2;
auto fal = fa + lev * nfc * nlat;
memset(fal, 0, nfc * nlat * sizeof(double));
for (long jm = 0; jm <= nt; jm++)
......@@ -413,8 +413,8 @@ sp2fctest(const double *sa, double *fa, const double *poli, long nlev, long nlat
const double sai = *sal++;
const double saris = sar * is;
const double saiis = sai * is;
double *restrict far = fal;
double *restrict fai = fal + nlat;
auto far = fal;
auto fai = fal + nlat;
#ifdef HAVE_OPENMP4
#pragma omp simd
#endif
......@@ -443,9 +443,9 @@ sp2fc(const double *sa, double *fa, const double *poli, long nlev, long nlat, lo
#endif
for (long lev = 0; lev < nlev; lev++)
{
const double *restrict pol = poli;
const double *restrict sal = sa + lev * nsp2;
double *fal = fa + lev * nfc * nlat;
auto pol = poli;
auto sal = sa + lev * nsp2;
auto fal = fa + lev * nfc * nlat;
memset(fal, 0, nfc * nlat * sizeof(double));
for (long jmm = 0; jmm <= nt; jmm++)
......@@ -474,7 +474,7 @@ sp2fc(const double *sa, double *fa, const double *poli, long nlev, long nlat, lo
}
void
fc2sp(double *fa, double *sa, const double *poli, long nlev, long nlat, long nfc, long nt)
fc2sp(const double *fa, double *sa, const double *poli, long nlev, long nlat, long nfc, long nt)
{
long nsp2 = (nt + 1) * (nt + 2);
......@@ -483,16 +483,16 @@ fc2sp(double *fa, double *sa, const double *poli, long nlev, long nlat, long nfc
#endif
for (long lev = 0; lev < nlev; lev++)
{
const double *restrict pol = poli;
double *fal = fa + lev * nfc * nlat;
double *sal = sa + lev * nsp2;
auto pol = poli;
auto fal = fa + lev * nfc * nlat;
auto sal = sa + lev * nsp2;
for (long jmm = 0; jmm <= nt; jmm++)
{
for (long jfc = jmm; jfc <= nt; jfc++)
{
const double *far = fal;
const double *fai = fal + nlat;
auto far = fal;
auto fai = fal + nlat;
double sar = 0.0;
double sai = 0.0;
#ifdef HAVE_OPENMP4
......
......@@ -22,8 +22,9 @@
#include <mpim_grid.h>
void
grid2spec(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
grid2spec(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
const auto &fcTrans = spTrans.fcTrans;
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDout);
const long nlon = gridInqXsize(gridIDin);
......@@ -33,17 +34,18 @@ grid2spec(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
std::vector<double> fpwork(nlat * nfc * nlev);
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
gp2fc(arrayIn, fpwork.data(), nlat, nlon, nfc);
else
gp2fc(sptrans.vtrig.data(), sptrans.ifax, arrayIn, fpwork.data(), nlat, nlon, nlev, nfc);
gp2fc(fcTrans.vtrig.data(), fcTrans.ifax, arrayIn, fpwork.data(), nlat, nlon, nlev, nfc);
fc2sp(fpwork.data(), arrayOut, sptrans.vpold.data(), nlev, nlat, nfc, ntr);
fc2sp(fpwork.data(), arrayOut, spTrans.vpold.data(), nlev, nlat, nfc, ntr);
}
void
spec2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
spec2grid(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
const auto &fcTrans = spTrans.fcTrans;
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDin);
const long nlon = gridInqXsize(gridIDout);
......@@ -53,29 +55,29 @@ spec2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
std::vector<double> fpwork(nlat * nfc * nlev);
sp2fc(arrayIn, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
sp2fc(arrayIn, fpwork.data(), spTrans.vpoli.data(), nlev, nlat, nfc, ntr);
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
fc2gp(fpwork.data(), arrayOut, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), arrayOut, nlat, nlon, nlev, nfc);
fc2gp(fcTrans.vtrig.data(), fcTrans.ifax, fpwork.data(), arrayOut, nlat, nlon, nlev, nfc);
}
void
four2spec(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
four2spec(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
(void) gridIDin;
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDout);
const long nlat = sptrans.nlat;
const long nlat = spTrans.nlat;
const long waves = ntr + 1;
const long nfc = waves * 2;
fc2sp(arrayIn, arrayOut, sptrans.vpold.data(), nlev, nlat, nfc, ntr);
fc2sp(arrayIn, arrayOut, spTrans.vpold.data(), nlev, nlat, nfc, ntr);
}
void
spec2four(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
spec2four(const SP_Transformation &spTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDin);
......@@ -84,11 +86,11 @@ spec2four(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
const long waves = ntr + 1;
nfc = waves * 2;
sp2fc(arrayIn, arrayOut, sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
sp2fc(arrayIn, arrayOut, spTrans.vpoli.data(), nlev, nlat, nfc, ntr);
}
void
four2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
four2grid(const FC_Transformation &fcTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDin);
......@@ -97,14 +99,14 @@ four2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
const long waves = ntr + 1;
const long nfc = waves * 2;
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
fc2gp(arrayIn, arrayOut, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
fc2gp(fcTrans.vtrig.data(), fcTrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
}
void
grid2four(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
grid2four(const FC_Transformation &fcTrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
const long nlev = 1;
const long ntr = gridInqTrunc(gridIDout);
......@@ -113,10 +115,10 @@ grid2four(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
const long waves = ntr + 1;
const long nfc = waves * 2;
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
gp2fc(arrayIn, arrayOut, nlat, nlon, nfc);
else
gp2fc(sptrans.vtrig.data(), sptrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
gp2fc(fcTrans.vtrig.data(), fcTrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
}
void
......@@ -137,7 +139,7 @@ speccut(int gridIDin, double *arrayIn, double *arrayOut, int *waves)
}
void
trans_uv2dv(SPTRANS &sptrans, long nlev, int gridID1, double *gu, double *gv, int gridID2, double *sd, double *svo)
trans_uv2dv(const SP_Transformation &spTrans, long nlev, int gridID1, double *gu, double *gv, int gridID2, double *sd, double *svo)
{
if (gridInqType(gridID1) != GRID_GAUSSIAN)
cdoAbort("unexpected grid1 type: %s instead of Gaussian", gridNamePtr(gridInqType(gridID1)));
......@@ -145,6 +147,7 @@ trans_uv2dv(SPTRANS &sptrans, long nlev, int gridID1, double *gu, double *gv, in
if (gridInqType(gridID2) != GRID_SPECTRAL)
cdoAbort("unexpected grid2 type: %s instead of spectral", gridNamePtr(gridInqType(gridID2)));
const auto &fcTrans = spTrans.fcTrans;
const long ntr = gridInqTrunc(gridID2);
const long nlon = gridInqXsize(gridID1);
const long nlat = gridInqYsize(gridID1);
......@@ -154,25 +157,25 @@ trans_uv2dv(SPTRANS &sptrans, long nlev, int gridID1, double *gu, double *gv, in
std::vector<double> fpwork1(nlat * nfc * nlev);
std::vector<double> fpwork2(nlat * nfc * nlev);
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
{
gp2fc(gu, fpwork1.data(), nlat, nlon, nfc);
gp2fc(gv, fpwork2.data(), nlat, nlon, nfc);
}
else
{
gp2fc(sptrans.vtrig.data(), sptrans.ifax, gu, fpwork1.data(), nlat, nlon, nlev, nfc);
gp2fc(sptrans.vtrig.data(), sptrans.ifax, gv, fpwork2.data(), nlat, nlon, nlev, nfc);
gp2fc(fcTrans.vtrig.data(), fcTrans.ifax, gu, fpwork1.data(), nlat, nlon, nlev, nfc);
gp2fc(fcTrans.vtrig.data(), fcTrans.ifax, gv, fpwork2.data(), nlat, nlon, nlev, nfc);
}
scaluv(fpwork1.data(), sptrans.vcoslat.data(), nlat, nfc * nlev);
scaluv(fpwork2.data(), sptrans.vcoslat.data(), nlat, nfc * nlev);
scaluv(fpwork1.data(), spTrans.vcoslat.data(), nlat, nfc * nlev);
scaluv(fpwork2.data(), spTrans.vcoslat.data(), nlat, nfc * nlev);
uv2dv(fpwork1.data(), fpwork2.data(), sd, svo, sptrans.vpol2.data(), sptrans.vpol3.data(), nlev, nlat, ntr);
uv2dv(fpwork1.data(), fpwork2.data(), sd, svo, spTrans.vpol2.data(), spTrans.vpol3.data(), nlev, nlat, ntr);
}
void
trans_dv2uv(SPTRANS &sptrans, DVTRANS &dvtrans, long nlev, int gridID1, double *sd, double *svo, int gridID2, double *gu,
trans_dv2uv(const SP_Transformation &spTrans, const DV_Transformation &dvTrans, long nlev, int gridID1, double *sd, double *svo, int gridID2, double *gu,
double *gv)
{
if (gridInqType(gridID1) != GRID_SPECTRAL)
......@@ -180,6 +183,7 @@ trans_dv2uv(SPTRANS &sptrans, DVTRANS &dvtrans, long nlev, int gridID1, double *
if (gridInqType(gridID2) != GRID_GAUSSIAN)
cdoWarning("unexpected grid2 type: %s instead of Gaussian", gridNamePtr(gridInqType(gridID2)));
const auto &fcTrans = spTrans.fcTrans;
const long ntr = gridInqTrunc(gridID1);
const long nlon = gridInqXsize(gridID2);
const long nlat = gridInqYsize(gridID2);
......@@ -190,23 +194,23 @@ trans_dv2uv(SPTRANS &sptrans, DVTRANS &dvtrans, long nlev, int gridID1, double *
double *su = gu;
double *sv = gv;
dv2uv(sd, svo, su, sv, dvtrans.f1.data(), dvtrans.f2.data(), ntr, dimsp, nlev);
dv2uv(sd, svo, su, sv, dvTrans.f1.data(), dvTrans.f2.data(), ntr, dimsp, nlev);
std::vector<double> fpwork(nlat * nfc * nlev);
sp2fc(su, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), sptrans.vrcoslat.data(), nlat, nfc * nlev);
sp2fc(su, fpwork.data(), spTrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), spTrans.vrcoslat.data(), nlat, nfc * nlev);
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
fc2gp(fpwork.data(), gu, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gu, nlat, nlon, nlev, nfc);
fc2gp(fcTrans.vtrig.data(), fcTrans.ifax, fpwork.data(), gu, nlat, nlon, nlev, nfc);
sp2fc(sv, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), sptrans.vrcoslat.data(), nlat, nfc * nlev);
sp2fc(sv, fpwork.data(), spTrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), spTrans.vrcoslat.data(), nlat, nfc * nlev);
if (sptrans.use_fftw)
if (fcTrans.use_fftw)
fc2gp(fpwork.data(), gv, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gv, nlat, nlon, nlev, nfc);
fc2gp(fcTrans.vtrig.data(), fcTrans.ifax, fpwork.data(), gv, nlat, nlon, nlev, nfc);
}
......@@ -17,19 +17,67 @@
#ifndef SPECSPACE_H
#define SPECSPACE_H
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <vector>
#include "cdo_options.h"
#include "transform.h"
class SPTRANS
class FC_Transformation
{
public:
bool use_fftw = false;
long nlon = 0;
long nlat = 0;
long ntr = 0;
long poldim = 0;
long ifax[10];
std::vector<double> vtrig;
FC_Transformation()
{
if (Options::Use_FFTW)