Commit 70fde1be authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

spectral transformation: changed address space from 32 to 64-bit integer.

parent 8c0afb5a
......@@ -3,6 +3,10 @@
* Using CDI library version 1.9.6
* Version 1.9.6 release
2018-10-16 Uwe Schulzweida
* spectral transformation: changed address space from 32 to 64-bit integer
2018-10-15 Uwe Schulzweida
* masklonlatbox: wrong result if lon1 > first lon || lon2 < last lon (bug introduce in 1.9.4) [Bug #8695]
......
......@@ -38,58 +38,54 @@
void *
Spectral(void *process)
{
int nrecs, nvars;
int nrecs;
int varID, levelID;
int index;
int gridIDsp = -1, gridIDgp = -1;
int gridID1 = -1, gridID2 = -1;
int gridID;
size_t nmiss;
int ncut = 0;
std::vector<int> waves;
int nlon, nlat, ntr;
SPTRANS *sptrans = NULL;
ListArray<int> listArrayInt;
cdoInitialize(process);
bool lcopy = UNCHANGED_RECORD;
const bool lcopy = UNCHANGED_RECORD;
// clang-format off
int GP2SP = cdoOperatorAdd("gp2sp", 0, 0, NULL);
int GP2SPL = cdoOperatorAdd("gp2spl", 0, 0, NULL);
int SP2GP = cdoOperatorAdd("sp2gp", 0, 0, NULL);
int SP2GPL = cdoOperatorAdd("sp2gpl", 0, 0, NULL);
int SP2SP = cdoOperatorAdd("sp2sp", 0, 0, NULL);
int SPCUT = cdoOperatorAdd("spcut", 0, 0, NULL);
const int GP2SP = cdoOperatorAdd("gp2sp", 0, 0, NULL);
const int GP2SPL = cdoOperatorAdd("gp2spl", 0, 0, NULL);
const int SP2GP = cdoOperatorAdd("sp2gp", 0, 0, NULL);
const int SP2GPL = cdoOperatorAdd("sp2gpl", 0, 0, NULL);
const int SP2SP = cdoOperatorAdd("sp2sp", 0, 0, NULL);
const int SPCUT = cdoOperatorAdd("spcut", 0, 0, NULL);
// clang-format on
int operatorID = cdoOperatorID();
const int operatorID = cdoOperatorID();
int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
const int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
int vlistID1 = cdoStreamInqVlist(streamID1);
int vlistID2 = vlistDuplicate(vlistID1);
const int vlistID1 = cdoStreamInqVlist(streamID1);
const int vlistID2 = vlistDuplicate(vlistID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
const int taxisID1 = vlistInqTaxis(vlistID1);
const int taxisID2 = taxisDuplicate(taxisID1);
vlistDefTaxis(vlistID2, taxisID2);
int ngrids = vlistNgrids(vlistID1);
/* find first spectral grid */
for (index = 0; index < ngrids; index++)
const int ngrids = vlistNgrids(vlistID1);
// find first spectral grid
for (int index = 0; index < ngrids; index++)
{
gridID = vlistGrid(vlistID1, index);
int gridID = vlistGrid(vlistID1, index);
if (gridInqType(gridID) == GRID_SPECTRAL)
{
gridIDsp = gridID;
break;
}
}
/* find first gaussian grid */
for (index = 0; index < ngrids; index++)
// find first gaussian grid
for (int index = 0; index < ngrids; index++)
{
gridID = vlistGrid(vlistID1, index);
int gridID = vlistGrid(vlistID1, index);
if (gridInqType(gridID) == GRID_GAUSSIAN)
{
gridIDgp = gridID;
......@@ -97,7 +93,7 @@ Spectral(void *process)
}
}
/* define output grid */
// define output grid
if (operatorID == GP2SP || operatorID == GP2SPL)
{
if (gridIDgp == -1) cdoWarning("No data on Gaussian grid found!");
......@@ -106,10 +102,10 @@ Spectral(void *process)
if (gridID1 != -1)
{
if (operatorID == GP2SP)
ntr = nlat_to_ntr(gridInqYsize(gridID1));
else
ntr = nlat_to_ntr_linear(gridInqYsize(gridID1));
const long nlon = gridInqXsize(gridID1);
const long nlat = gridInqYsize(gridID1);
long ntr = (operatorID == GP2SP) ? nlat_to_ntr(nlat) : nlat_to_ntr_linear(nlat);
if (gridIDsp != -1)
if (ntr != gridInqTrunc(gridIDsp)) gridIDsp = -1;
......@@ -128,10 +124,7 @@ Spectral(void *process)
gridID2 = gridIDsp;
nlon = gridInqXsize(gridID1);
nlat = gridInqYsize(gridID1);
ntr = gridInqTrunc(gridID2);
sptrans = sptrans_new(nlon, nlat, ntr, 0);
}
}
......@@ -145,31 +138,23 @@ Spectral(void *process)
{
if (gridIDgp != -1)
{
if (operatorID == SP2GP)
ntr = nlat_to_ntr(gridInqYsize(gridIDgp));
else
ntr = nlat_to_ntr_linear(gridInqYsize(gridIDgp));
const long nlat = gridInqYsize(gridIDgp);
const long ntr = (operatorID == SP2GP) ? nlat_to_ntr(nlat) : nlat_to_ntr_linear(nlat);
if (gridInqTrunc(gridIDsp) != ntr) gridIDgp = -1;
}
if (gridIDgp == -1)
{
char gridname[20];
if (operatorID == SP2GP)
snprintf(gridname, sizeof(gridname), "t%dgrid", gridInqTrunc(gridIDsp));
else
snprintf(gridname, sizeof(gridname), "tl%dgrid", gridInqTrunc(gridIDsp));
snprintf(gridname, sizeof(gridname), "t%s%dgrid", (operatorID == SP2GP)?"":"l", gridInqTrunc(gridIDsp));
gridIDgp = grid_from_name(gridname);
}
gridID2 = gridIDgp;
ntr = gridInqTrunc(gridID1);
nlon = gridInqXsize(gridID2);
nlat = gridInqYsize(gridID2);
const long ntr = gridInqTrunc(gridID1);
const long nlon = gridInqXsize(gridID2);
const long nlat = gridInqYsize(gridID2);
sptrans = sptrans_new(nlon, nlat, ntr, 0);
}
}
......@@ -181,8 +166,8 @@ Spectral(void *process)
if (gridID1 != -1)
{
if (!isdigit(operatorArgv()[0][0])) cdoAbort("parameter truncation must comprise only digits [0-9]!");
int ntr = parameter2int(operatorArgv()[0]);
int nsp = (ntr + 1) * (ntr + 2);
const long ntr = parameter2int(operatorArgv()[0]);
const long nsp = (ntr + 1) * (ntr + 2);
gridIDsp = gridCreate(GRID_SPECTRAL, nsp);
gridDefTrunc(gridIDsp, ntr);
gridDefComplexPacking(gridIDsp, 1);
......@@ -194,20 +179,19 @@ Spectral(void *process)
}
else if (operatorID == SPCUT)
{
long i, j, maxntr;
gridID1 = gridIDsp;
operatorInputArg("wave numbers");
if (gridID1 != -1)
{
maxntr = 1 + gridInqTrunc(gridID1);
ncut = listArrayInt.argvToInt(operatorArgc(), operatorArgv());
const long maxntr = 1 + gridInqTrunc(gridID1);
const long ncut = listArrayInt.argvToInt(operatorArgc(), operatorArgv());
const int *wnums = listArrayInt.data();
waves.resize(maxntr);
for (i = 0; i < maxntr; i++) waves[i] = 1;
for (i = 0; i < ncut; i++)
for (long i = 0; i < maxntr; i++) waves[i] = 1;
for (long i = 0; i < ncut; i++)
{
j = wnums[i] - 1;
const long j = wnums[i] - 1;
if (j < 0 || j >= maxntr) cdoAbort("wave number %ld out of range (min=1, max=%l qd)!", wnums[i], maxntr);
waves[j] = 0;
}
......@@ -218,7 +202,7 @@ Spectral(void *process)
gridID2 = gridIDsp;
}
nvars = vlistNvars(vlistID2);
const int nvars = vlistNvars(vlistID2);
std::vector<bool> vars(nvars);
for (varID = 0; varID < nvars; varID++)
{
......@@ -227,10 +211,10 @@ Spectral(void *process)
if (gridID1 != -1) vlistChangeGrid(vlistID2, gridID1, gridID2);
int streamID2 = cdoStreamOpenWrite(cdoStreamName(1), cdoFiletype());
const int streamID2 = cdoStreamOpenWrite(cdoStreamName(1), cdoFiletype());
pstreamDefVlist(streamID2, vlistID2);
size_t gridsizemax = vlistGridsizeMax(vlistID1);
const size_t gridsizemax = vlistGridsizeMax(vlistID1);
std::vector<double> array1(gridsizemax);
std::vector<double> array2;
if (gridID2 != -1) array2.resize(gridInqSize(gridID2));
......
......@@ -56,14 +56,14 @@ dv2ps(const double *restrict div, double *restrict pot, long nlev, long ntr)
}
void
dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, int nt, int nsp, int nlev)
dv2uv(double *d, double *o, double *u, double *v, double *f, 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 */
/* f(nsp/2) , g(nsp/2) ! factor tables */
int l, m, n;
int i;
long l, m, n;
long i;
for (l = 0; l < nlev; l++)
{
......@@ -160,9 +160,9 @@ dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, int nt,
}
/*
void scaluv(double *fu, double rclat[], int nlat, int lot)
void scaluv(double *fu, double rclat[], long nlat, long lot)
{
int l, lat;
long l, lat;
double *ful;
#if defined (_OPENMP)
......@@ -180,12 +180,10 @@ void scaluv(double *fu, double rclat[], int nlat, int lot)
*/
void
scaluv(double *fu, double *rclat, int nlat, int lot)
scaluv(double *fu, double *rclat, long nlat, long lot)
{
int l, lat;
for (l = 0; l < lot; l++)
for (lat = 0; lat < nlat; lat++)
for (long l = 0; l < lot; l++)
for (long lat = 0; lat < nlat; lat++)
{
*fu *= rclat[lat];
fu++;
......@@ -193,16 +191,16 @@ scaluv(double *fu, double *rclat, int nlat, int lot)
}
void
uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, int klev, int nlat, int nt)
uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, long klev, long nlat, long nt)
{
int lev, jmm, jfc, lat, nfc, nsp2;
long lev, jmm, jfc, lat;
double dir, dii, vor, voi;
double *ufr, *ufi, *vfr, *vfi;
double *ful, *fvl, *sdl, *svl;
double *po2, *po3;
nsp2 = (nt + 1) * (nt + 2);
nfc = (nt + 1) * 2;
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, \
......
......@@ -38,7 +38,7 @@ void gaussaw(double *pa, double *pw, size_t nlat);
}
static void
jspleg1(double *pleg, double plat, int ktrunc, double *work)
jspleg1(double *pleg, double plat, long ktrunc, double *work)
{
/*
jspleg1 - Routine to calculate legendre functions
......@@ -89,7 +89,7 @@ jspleg1(double *pleg, double plat, int ktrunc, double *work)
None
*/
int itout1, i1m, ilm, jm, jcn, im2;
long itout1, i1m, ilm, jm, jcn, im2;
double zsin, zcos, zf1m, zre1, zf2m, znsqr, ze1, ze2;
double zjmsqr;
double *zhlp1, *zhlp2, *zhlp3;
......@@ -173,11 +173,11 @@ jspleg1(double *pleg, double plat, int ktrunc, double *work)
/* and their meridional derivatives */
/* ============================================= */
static void
phcs(double *pnm, double *hnm, int waves, double pmu, double *ztemp1, double *ztemp2)
phcs(double *pnm, double *hnm, long waves, double pmu, double *ztemp1, double *ztemp2)
{
int twowaves;
long twowaves;
int jk, jn, jm;
long jk, jn, jm;
double jnmjk;
double zcos2;
......@@ -287,15 +287,15 @@ phcs(double *pnm, double *hnm, int waves, double pmu, double *ztemp1, double *zt
}
void
after_legini_full(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict pdev, double *restrict pol2,
after_legini_full(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict pdev, double *restrict pol2,
double *restrict pol3, double *restrict coslat)
{
int jgl, jm, jn;
int jsp;
long jgl, jm, jn;
long jsp;
double gmusq;
int waves = ntr + 1;
int dimsp = (ntr + 1) * (ntr + 2);
long waves = ntr + 1;
long dimsp = (ntr + 1) * (ntr + 2);
std::vector<double> gmu(nlat);
std::vector<double> gwt(nlat);
......@@ -349,12 +349,12 @@ after_legini_full(int ntr, int nlat, double *restrict poli, double *restrict pol
}
void
after_legini(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict coslat)
after_legini(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict coslat)
{
int is, lats;
long is, lats;
int waves = ntr + 1;
int dimpnm = (ntr + 1) * (ntr + 4) / 2;
long waves = ntr + 1;
long dimpnm = (ntr + 1) * (ntr + 4) / 2;
std::vector<double> gmu(nlat);
std::vector<double> gwt(nlat);
......@@ -362,19 +362,19 @@ after_legini(int ntr, int nlat, double *restrict poli, double *restrict pold, do
std::vector<double> work(3 * waves);
gaussaw(gmu.data(), gwt.data(), nlat);
for (int jgl = 0; jgl < nlat; ++jgl) gwt[jgl] *= 0.5;
for (long jgl = 0; jgl < nlat; ++jgl) gwt[jgl] *= 0.5;
for (int jgl = 0; jgl < nlat; ++jgl) coslat[jgl] = sqrt(1.0 - gmu[jgl] * gmu[jgl]);
for (long jgl = 0; jgl < nlat; ++jgl) coslat[jgl] = sqrt(1.0 - gmu[jgl] * gmu[jgl]);
for (int jgl = 0; jgl < nlat / 2; jgl++)
for (long jgl = 0; jgl < nlat / 2; jgl++)
{
double zgwt = gwt[jgl];
jspleg1(pnm.data(), gmu[jgl], ntr, work.data());
int latn = jgl;
int isp = 0;
for (int jm = 0; jm < waves; jm++)
long latn = jgl;
long isp = 0;
for (long jm = 0; jm < waves; jm++)
{
#if defined(SX)
#pragma vdir nodep
......@@ -382,7 +382,7 @@ after_legini(int ntr, int nlat, double *restrict poli, double *restrict pold, do
#ifdef HAVE_OPENMP4
#pragma omp simd
#endif
for (int jn = 0; jn < waves - jm; jn++)
for (long jn = 0; jn < waves - jm; jn++)
{
is = (jn + 1) % 2 * 2 - 1;
lats = latn - jgl + nlat - jgl - 1;
......@@ -492,14 +492,14 @@ sp2fc(const double *sa, double *fa, const double *poli, long nlev, long nlat, lo
}
void
fc2sp(double *fa, double *sa, double *poli, int nlev, int nlat, int nfc, int nt)
fc2sp(double *fa, double *sa, double *poli, long nlev, long nlat, long nfc, long nt)
{
int nsp2 = (nt + 1) * (nt + 2);
long nsp2 = (nt + 1) * (nt + 2);
#ifdef _OPENMP
#pragma omp parallel for default(shared)
#endif
for (int lev = 0; lev < nlev; lev++)
for (long lev = 0; lev < nlev; lev++)
{
const double *restrict pol = poli;
double *fal = fa + lev * nfc * nlat;
......@@ -508,7 +508,7 @@ fc2sp(double *fa, double *sa, double *poli, int nlev, int nlat, int nfc, int nt)
const double *restrict far;
const double *restrict fai;
double sar, sai;
int jmm, jfc, lat;
long jmm, jfc, lat;
for (jmm = 0; jmm <= nt; jmm++)
{
......@@ -540,9 +540,9 @@ fc2sp(double *fa, double *sa, double *poli, int nlev, int nlat, int nfc, int nt)
/* ======================================== */
void
sp2sp(double *arrayIn, int truncIn, double *arrayOut, int truncOut)
sp2sp(double *arrayIn, long truncIn, double *arrayOut, long truncOut)
{
int n, m;
long n, m;
if (truncOut <= truncIn)
{
......@@ -585,9 +585,9 @@ sp2sp(double *arrayIn, int truncIn, double *arrayOut, int truncOut)
/* ======================================== */
void
spcut(double *arrayIn, double *arrayOut, int trunc, int *waves)
spcut(double *arrayIn, double *arrayOut, long trunc, int *waves)
{
int n, m;
long n, m;
for (n = 0; n <= trunc; n++)
{
......
......@@ -109,9 +109,9 @@ struct Control
int Truncation;
int Waves;
int Dim3FC, Dim3SP, Dim3GP;
int DimFC, DimGP, DimSP;
int DimSP_half;
long Dim3FC, Dim3SP, Dim3GP;
long DimFC, DimGP, DimSP;
long DimSP_half;
double *poli;
double *pold;
......@@ -176,15 +176,15 @@ void fc2gp(double * trig, long * ifax, double * fc, double * gp, long nlat, long
void gp2fc(double *trig, long *ifax, const double * gp, double * fc, long nlat, long nlon, long nlev, long nfc);
/* Convert Spectral Array to new resolution */
void sp2sp(double *arrayIn, int truncIn, double *arrayOut, int truncOut);
void sp2sp(double *arrayIn, long truncIn, double *arrayOut, long truncOut);
void sp2fc(const double *sa, double *fa, const double *poli, long nlev, long nlat, long nfc, long nt);
void fc2sp(double *fa, double *sa, double *poli, int klev, int nlat, int nfc, int nt);
void fc2sp(double *fa, double *sa, double *poli, long klev, long nlat, long nfc, long nt);
/* Physc */
void dv2ps(const double * div, double * pot, long nlev, long ntr);
void dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, int nt, int nsp, int nlev);
void scaluv(double *fu, double rclat[], int nlat, int lot);
void uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, int klev, int nlat, int nt);
void dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, long nt, long nsp, long nlev);
void scaluv(double *fu, 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);
void geninx(long ntr, double *f, double *g);
/* clang-format off */
......
......@@ -2717,17 +2717,17 @@ after_EchamDependencies(struct Variable *vars, int ncodes, int type, int source)
CheckDependencies(vars, 0);
}
void after_legini_full(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict pdev,
void after_legini_full(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict pdev,
double *restrict pol2, double *restrict pol3, double *restrict coslat);
void after_legini(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict coslat);
void after_legini(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict coslat);
void
after_legini_setup(struct Control *globs, struct Variable *vars)
{
int ntr = globs->Truncation;
int nlat = globs->Latitudes;
int dimsp = (ntr + 1) * (ntr + 2);
int pdim = dimsp / 2 * nlat;
long ntr = globs->Truncation;
long nlat = globs->Latitudes;
long dimsp = (ntr + 1) * (ntr + 2);
long pdim = dimsp / 2 * nlat;
globs->poli = (double *) Malloc(pdim * sizeof(double));
......@@ -2772,7 +2772,7 @@ after_legini_setup(struct Control *globs, struct Variable *vars)
else
after_legini(ntr, nlat, globs->poli, globs->pold, globs->coslat);
for (int jgl = 0; jgl < nlat; ++jgl) globs->rcoslat[jgl] = 1.0 / globs->coslat[jgl];
for (long jgl = 0; jgl < nlat; ++jgl) globs->rcoslat[jgl] = 1.0 / globs->coslat[jgl];
for (int jgl = 0; jgl < nlat; ++jgl) globs->DerivationFactor[jgl] = globs->rcoslat[jgl] / PlanetRadius;
for (long jgl = 0; jgl < nlat; ++jgl) globs->DerivationFactor[jgl] = globs->rcoslat[jgl] / PlanetRadius;
}
......@@ -293,6 +293,25 @@ arrayNumMV(const size_t len, const double *restrict array, const double missval)
return nmiss;
}
size_t
arrayNumMV_f(const size_t len, const float *restrict array, const double missval)
{
size_t nmiss = 0;
if (DBL_IS_NAN(missval))
{
for (size_t i = 0; i < len; ++i)
if (DBL_IS_EQUAL(array[i], missval)) nmiss++;
}
else
{
for (size_t i = 0; i < len; ++i)
if (IS_EQUAL(array[i], missval)) nmiss++;
}
return nmiss;
}
double
arrayMin(const size_t len, const double *restrict array)
{
......
......@@ -109,6 +109,7 @@ arrayCopy(const size_t len, const T *array1, T * array2)
}
size_t arrayNumMV(const size_t len, const double *array, const double missval);
size_t arrayNumMV_f(const size_t len, const float *array, const double missval);
double arrayMin(const size_t len, const double *array);
double arrayMax(const size_t len, const double *array);
......
......@@ -21,22 +21,22 @@
void geninx(long ntr, double *f, double *g);
void scaluv(double *fu, double *rclat, int nlat, int lot);
void uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, int klev, int nlat, int nt);
void dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, int nt, int nsp, int nlev);
void uv2dv(double *fu, double *fv, double *sd, double *sv, double *pol2, double *pol3, long klev, long nlat, long nt);
void dv2uv(double *d, double *o, double *u, double *v, double *f, double *g, long nt, long nsp, long nlev);
void after_legini_full(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict pdev,
void after_legini_full(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict pdev,
double *restrict pol2, double *restrict pol3, double *restrict coslat);
void after_legini(int ntr, int nlat, double *restrict poli, double *restrict pold, double *restrict coslat);
void after_legini(long ntr, long nlat, double *restrict poli, double *restrict pold, double *restrict coslat);
void
grid2spec(SPTRANS *sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
int nlev = 1;
int ntr = gridInqTrunc(gridIDout);
int nlon = gridInqXsize(gridIDin);
int nlat = gridInqYsize(gridIDin);
int waves = ntr + 1;
int nfc = waves * 2;
long nlev = 1;
long ntr = gridInqTrunc(gridIDout);
long nlon = gridInqXsize(gridIDin);
long nlat = gridInqYsize(gridIDin);
long waves = ntr + 1;
long nfc = waves * 2;
std::vector<double> fpwork(nlat * nfc * nlev);
......@@ -47,12 +47,12 @@ grid2spec(SPTRANS *sptrans, int gridIDin, double *arrayIn, int gridIDout, double
void
spec2grid(SPTRANS *sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
int nlev = 1;
int ntr = gridInqTrunc(gridIDin);
int nlon = gridInqXsize(gridIDout);
int nlat = gridInqYsize(gridIDout);
int waves = ntr + 1;
int nfc = waves * 2;
long nlev = 1;
long ntr = gridInqTrunc(gridIDin);
long nlon = gridInqXsize(gridIDout);
long nlat = gridInqYsize(gridIDout);
long waves = ntr + 1;
long nfc = waves * 2;
std::vector<double> fpwork(nlat * nfc * nlev);
......@@ -64,11 +64,11 @@ void
four2spec(SPTRANS *sptrans, int gridIDin, double *arrayIn, int gridIDout, double *arrayOut)
{
(void) gridIDin;