Commit a95a47ac authored by Oliver Heidmann's avatar Oliver Heidmann
Browse files

Merge branch 'develop' of git.mpimet.mpg.de:cdo into develop

parents 4d33e781 233d9cad
2017-08-14 Uwe Schulzweida
* eof, eof3d: set default of env. CDO_WEIGHT_MODE to off
* eof3d: preserve variable name on first output file [report: Frank Kauker]
2017-08-12 Uwe Schulzweida
* eof3d: weight was allocated for only one level (bug fix) [report: Frank Kauker]
2017-07-27 Uwe Schulzweida
* Using CDI library version 1.9.0
......
CDO NEWS
--------
Version 1.9.1 (27 September 2017):
Changes operators:
* eof, eof3d: set default of environment variable CDO_WEIGHT_MODE to off
Fixed bugs:
* eof3d: weight array was allocated for only one level (bug fix)
Version 1.9.0 (27 July 2017):
New features:
......
......@@ -198,7 +198,7 @@ keepaspectratio]{cdo_libdep.pdf}}%
\end{picture}
\begin{flushright}
\large \textbf{Climate Data Operators \\ Version 1.9.0 \\ August 2017}
\large \textbf{Climate Data Operators \\ Version 1.9.1 \\ September 2017}
\end{flushright}
\vfill
......
......@@ -14,7 +14,7 @@
\put(0,0.0){\line(1,0){3.95}}
\end{picture}
\begin{flushright}
{\small{Climate Data Operators \\ Version 1.9.0 \\ August 2017}}
{\small{Climate Data Operators \\ Version 1.9.1 \\ September 2017}}
\end{flushright}
\vspace*{0mm}
......
......@@ -73,7 +73,7 @@ Is used to choose the algorithm for eigenvalue calculation. Options are 'jacobi'
a one-sided parallel jacobi-algorithm (only executed in parallel if -P flag is set)
and 'danielson_lanczos' for a non-parallel d/l algorithm. The default setting is 'jacobi'.
@Item = CDO_WEIGHT_MODE
It is used to set the weight mode. The default is 'on'. Set it to 'off' for a non weighted version.
It is used to set the weight mode. The default is 'off'. Set it to 'on' for a weighted version.
@Item = MAX_JACOBI_ITER
Is the maximum integer number of annihilation sweeps that is executed if the
jacobi-algorithm is used to compute the eigen values. The default value is 12.
......
libcdi @ 405563ff
Subproject commit b6dd1e7488d2fe2c3d51f2858b73b848611d561b
Subproject commit 405563ffb3cd44fc1abc148bac4f306f2ebc4ca8
......@@ -43,20 +43,20 @@
// NO MISSING VALUE SUPPORT ADDED SO FAR
static
void scale_eigvec_grid(double *restrict out, int tsID, int npack, const int *restrict pack, const double *restrict weight, double **covar, double sum_w)
void scale_eigvec_grid(double *restrict out, int tsID, size_t npack, const size_t *restrict pack, const double *restrict weight, double **covar, double sum_w)
{
for ( int i = 0; i < npack; ++i )
for ( size_t i = 0; i < npack; ++i )
out[pack[i]] = covar[tsID][i] / sqrt(weight[pack[i]]/sum_w);
}
static
void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const int *restrict pack, const double *restrict weight,
void scale_eigvec_time(double *restrict out, int tsID, int nts, size_t npack, const size_t *restrict pack, const double *restrict weight,
double **covar, double **data, double missval, double sum_w)
{
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(npack, nts, tsID, pack, data, covar, out)
#endif
for ( int i = 0; i < npack; ++i )
for ( size_t i = 0; i < npack; ++i )
{
double sum = 0;
for ( int j = 0; j < nts; ++j )
......@@ -65,9 +65,9 @@ void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const
out[pack[i]] = sum;
}
/*
for ( int j = 0; j < nts; ++j )
for ( size_t j = 0; j < nts; ++j )
{
for ( int i = 0; i < npack; ++i )
for ( size_t i = 0; i < npack; ++i )
out[pack[i]] += data[j][i] * covar[tsID][j];
}
*/
......@@ -79,7 +79,7 @@ void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const
#pragma omp parallel for default(none) reduction(+:sum) \
shared(out,weight,pack,npack)
#endif
for ( int i = 0; i < npack; ++i )
for ( size_t i = 0; i < npack; ++i )
{
// do not need to account for weights as eigenvectors are non-weighted
sum += weight[pack[i]] * out[pack[i]] * out[pack[i]];
......@@ -91,14 +91,14 @@ void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(npack,pack,sum,out)
#endif
for ( int i = 0; i < npack; ++i ) out[pack[i]] /= sum;
for ( size_t i = 0; i < npack; ++i ) out[pack[i]] /= sum;
}
else
{
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(npack,pack,out,missval)
#endif
for ( int i = 0; i < npack; ++i ) out[pack[i]] = missval;
for ( size_t i = 0; i < npack; ++i ) out[pack[i]] = missval;
}
}
......@@ -141,7 +141,7 @@ enum T_EIGEN_MODE get_eigenmode(void)
enum T_WEIGHT_MODE get_weightmode(void)
{
enum T_WEIGHT_MODE weight_mode = WEIGHT_ON;
enum T_WEIGHT_MODE weight_mode = WEIGHT_OFF;
char *envstr = getenv("CDO_WEIGHT_MODE");
if ( envstr )
......@@ -171,7 +171,7 @@ void *EOFs(void * argument)
int nmiss;
int varID, levelID;
int nts = 0;
int n = 0;
size_t n = 0;
int grid_space = 0, time_space = 0;
int timer_cov = 0, timer_eig = 0;
......@@ -217,7 +217,7 @@ void *EOFs(void * argument)
int vlistID1 = pstreamInqVlist(streamID1);
int taxisID1 = vlistInqTaxis(vlistID1);
int gridID1 = vlistInqVarGrid(vlistID1, 0);
int gridsize = vlistGridsizeMax(vlistID1);
size_t gridsize = vlistGridsizeMax(vlistID1);
int nvars = vlistNvars(vlistID1);
int nrecs;
......@@ -226,19 +226,6 @@ void *EOFs(void * argument)
if ( vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index))
cdoAbort("Too many different grids!");
double *weight = (double *) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; ++i ) weight[i] = 1.;
if ( weight_mode == WEIGHT_ON )
{
int wstatus = gridWeights(gridID1, weight);
if ( wstatus != 0 )
{
weight_mode = WEIGHT_OFF;
cdoWarning("Using constant grid cell area weights!");
}
}
/* eigenvalues */
/* COUNT NUMBER OF TIMESTEPS if EOF_ or EOF_TIME */
......@@ -253,17 +240,17 @@ void *EOFs(void * argument)
while ( pstreamInqTimestep(streamID1, nts) ) nts++;
if ( cdoVerbose ) cdoPrint("Counted %i timeSteps", nts);
pstreamClose(streamID1);
streamID1 = pstreamOpenRead(cdoStreamName(0));
vlistID1 = pstreamInqVlist(streamID1);
taxisID1 = vlistInqTaxis(vlistID1);
}
else
if ( cdoVerbose ) cdoPrint("Found %i timeSteps", nts);
pstreamClose(streamID1);
streamID1 = pstreamOpenRead(cdoStreamName(0));
vlistID1 = pstreamInqVlist(streamID1);
taxisID1 = vlistInqTaxis(vlistID1);
if ( nts < gridsize || operfunc == EOF_TIME )
if ( (size_t)nts < gridsize || operfunc == EOF_TIME )
{
time_space = 1;
grid_space = 0;
......@@ -297,11 +284,11 @@ void *EOFs(void * argument)
{
if ( ((double)gridsize)*gridsize > (double)LONG_MAX ) cdoAbort("Grid space too large!");
if ( n_eig > gridsize )
if ( (size_t)n_eig > gridsize )
{
cdoWarning("Solving in spatial space");
cdoWarning("Number of eigen-functions to write out is bigger than grid size");
cdoWarning("Setting n_eig to %d", gridsize);
cdoWarning("Setting n_eig to %zu", gridsize);
cdoWarning("If You want to force a solution in time-space use operator eoftime");
n_eig = gridsize;
}
......@@ -312,9 +299,22 @@ void *EOFs(void * argument)
cdoPrint("Calculating %d eigenvectors and %d eigenvalues in %s",
n_eig, n, grid_space==1?"grid_space" : "time_space");
double *weight = (double *) Malloc(gridsize*sizeof(double));
for ( size_t i = 0; i < gridsize; ++i ) weight[i] = 1.;
if ( weight_mode == WEIGHT_ON )
{
int wstatus = gridWeights(gridID1, weight);
if ( wstatus != 0 )
{
weight_mode = WEIGHT_OFF;
cdoWarning("Using constant grid cell area weights!");
}
}
/* allocation of temporary fields and output structures */
int npack = -1;
int *pack = (int *) Malloc(gridsize*sizeof(int));
size_t npack = ULONG_MAX;
size_t *pack = (size_t *) Malloc(gridsize*sizeof(size_t));
double *in = (double *) Malloc(gridsize*sizeof(double));
eofdata_t **eofdata = (eofdata_t **) Malloc(nvars*sizeof(eofdata_t*));
......@@ -341,7 +341,7 @@ void *EOFs(void * argument)
}
if ( cdoVerbose )
cdoPrint("Allocated eigenvalue/eigenvector structures with nts=%d gridsize=%d", nts, gridsize);
cdoPrint("Allocated eigenvalue/eigenvector structures with nts=%d gridsize=%zu", nts, gridsize);
double *covar_array = NULL;
double **covar = NULL;
......@@ -362,10 +362,10 @@ void *EOFs(void * argument)
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
missval = vlistInqVarMissval(vlistID1, varID);
if ( npack == -1 )
if ( npack == ULONG_MAX )
{
npack = 0;
for ( int i = 0; i < gridsize; ++i )
for ( size_t i = 0; i < gridsize; ++i )
{
if ( !DBL_IS_EQUAL(weight[i], 0.0) && !DBL_IS_EQUAL(weight[i], missval) &&
!DBL_IS_EQUAL(in[i], missval) )
......@@ -378,12 +378,12 @@ void *EOFs(void * argument)
if ( weight_mode == WEIGHT_ON )
{
sum_w = 0;
for ( int i = 0; i < npack; i++ ) sum_w += weight[pack[i]];
for ( size_t i = 0; i < npack; i++ ) sum_w += weight[pack[i]];
}
}
int ipack = 0;
for ( int i = 0; i < gridsize; ++i )
size_t ipack = 0;
for ( size_t i = 0; i < gridsize; ++i )
{
if ( !DBL_IS_EQUAL(weight[i], 0.0) && !DBL_IS_EQUAL(weight[i], missval) &&
!DBL_IS_EQUAL(in[i], missval) )
......@@ -398,11 +398,11 @@ void *EOFs(void * argument)
if ( !eofdata[varID][levelID].init )
{
n = npack;
double *covar_array = (double *) Malloc(((size_t)npack)*npack*sizeof(double));
double *covar_array = (double *) Malloc(npack*npack*sizeof(double));
covar = (double **) Malloc(npack*sizeof(double *));
for ( int i = 0; i < npack; ++i ) covar[i] = covar_array + ((size_t)npack)*i;
for ( int i = 0; i < npack; ++i )
for ( int j = 0; j < npack; ++j ) covar[i][j] = 0;
for ( size_t i = 0; i < npack; ++i ) covar[i] = covar_array + npack*i;
for ( size_t i = 0; i < npack; ++i )
for ( size_t j = 0; j < npack; ++j ) covar[i][j] = 0;
eofdata[varID][levelID].covar_array = covar_array;
eofdata[varID][levelID].covar = covar;
......@@ -414,18 +414,16 @@ void *EOFs(void * argument)
#if defined(_OPENMP)
#pragma omp parallel for default(shared)
#endif
for ( int ipack = 0; ipack < npack; ++ipack )
{
for ( int jpack = ipack; jpack < npack; ++jpack )
covar[ipack][jpack] += in[pack[ipack]] * in[pack[jpack]];
}
for ( size_t ipack = 0; ipack < npack; ++ipack )
for ( size_t jpack = ipack; jpack < npack; ++jpack )
covar[ipack][jpack] += in[pack[ipack]] * in[pack[jpack]];
}
else if ( time_space )
{
double *data = (double *) Malloc(npack*sizeof(double));
eofdata[varID][levelID].data[tsID] = data;
for ( int ipack = 0; ipack < npack; ipack++ )
for ( size_t ipack = 0; ipack < npack; ipack++ )
data[ipack] = in[pack[ipack]];
}
......@@ -479,7 +477,7 @@ void *EOFs(void * argument)
double *eig_val = NULL;
int nts_out = nts;
if ( npack < nts ) nts_out = npack;
if ( npack < (size_t)nts ) nts_out = npack;
for ( tsID = 0; tsID < nts_out; tsID++ )
{
......@@ -529,11 +527,11 @@ void *EOFs(void * argument)
covar = eofdata[varID][levelID].covar;
for ( int ipack = 0; ipack < npack; ++ipack )
for ( size_t ipack = 0; ipack < npack; ++ipack )
{
int j;
int i = pack[ipack];
for ( int jpack = 0; jpack < npack; ++jpack )
size_t j;
size_t i = pack[ipack];
for ( size_t jpack = 0; jpack < npack; ++jpack )
{
if ( jpack < ipack )
{
......@@ -553,7 +551,7 @@ void *EOFs(void * argument)
else if ( time_space )
{
if ( cdoVerbose )
cdoPrint("allocating covar with %i x %i elements | npack=%i", nts, nts, npack);
cdoPrint("allocating covar with %i x %i elements | npack=%zu", nts, nts, npack);
covar_array = (double *) Malloc(nts*nts*sizeof(double));
covar = (double **) Malloc(nts*sizeof(double *));
......@@ -575,7 +573,7 @@ void *EOFs(void * argument)
{
double *df2p = data[j2];
double sum = 0;
for ( int i = 0; i < npack; i++ )
for ( size_t i = 0; i < npack; i++ )
sum += weight[pack[i]]*df1p[i]*df2p[i];
covar[j1][j2] = sum / sum_w / nts;
}
......@@ -599,7 +597,7 @@ void *EOFs(void * argument)
if ( cdoTimer ) timer_stop(timer_eig);
/* NOW: covar contains the eigenvectors, eig_val the eigenvalues */
for ( int i = 0; i < gridsize; ++i ) out[i] = missval;
for ( size_t i = 0; i < gridsize; ++i ) out[i] = missval;
// for ( int i = 0; i < n; i++ ) eig_val[i] *= sum_w;
} // first_call
......@@ -615,7 +613,7 @@ void *EOFs(void * argument)
else if ( time_space ) scale_eigvec_time(out, tsID, nts, npack, pack, weight, covar, data, missval, sum_w);
nmiss = 0;
for ( int i = 0; i < gridsize; i++ ) if ( DBL_IS_EQUAL(out[i], missval) ) nmiss++;
for ( size_t i = 0; i < gridsize; i++ ) if ( DBL_IS_EQUAL(out[i], missval) ) nmiss++;
pstreamDefRecord(streamID3, varID, levelID);
pstreamWriteRecord(streamID3, out, nmiss);
......
......@@ -53,8 +53,10 @@ void *EOF3d(void * argument)
size_t temp_size = 0, npack = 0;
int varID, levelID;
int missval_warning = 0;
int nmiss, ngrids, n = 0, nlevs = 0;
int offset;
int nmiss, ngrids;
int n = 0;
size_t nlevs = 0;
size_t offset;
int timer_cov = 0, timer_eig = 0;
int calendar = CALENDAR_STANDARD;
......@@ -88,31 +90,14 @@ void *EOF3d(void * argument)
enum T_EIGEN_MODE eigen_mode = get_eigenmode();
enum T_WEIGHT_MODE weight_mode = get_weightmode();
int streamID1 = pstreamOpenRead(cdoStreamName(0));
int vlistID1 = pstreamInqVlist(streamID1);
int gridID1 = vlistInqVarGrid(vlistID1, 0);
long gridsize = vlistGridsizeMax(vlistID1);
int nvars = vlistNvars(vlistID1);
int nrecs;
double *weight = (double *) Malloc(gridsize*sizeof(double));
for ( int i = 0; i < gridsize; ++i ) weight[i] = 1.;
if ( weight_mode == WEIGHT_ON )
{
int wstatus = gridWeights(gridID1, weight);
if ( wstatus != 0 )
{
weight_mode = WEIGHT_OFF;
cdoWarning("Using constant grid cell area weights!");
}
}
/* eigenvalues */
if ( operfunc == EOF3D_SPATIAL )
cdoAbort("Operator not Implemented - use eof3d or eof3dtime instead");
int streamID1 = pstreamOpenRead(cdoStreamName(0));
int vlistID1 = pstreamInqVlist(streamID1);
/* COUNT NUMBER OF TIMESTEPS if EOF3D_ or EOF3D_TIME */
int nts = vlistNtsteps(vlistID1);
if ( nts == -1 )
......@@ -121,14 +106,15 @@ void *EOF3d(void * argument)
while ( pstreamInqTimestep(streamID1, nts) ) nts++;
if ( cdoVerbose ) cdoPrint("Counted %i timeSteps", nts);
pstreamClose(streamID1);
streamID1 = pstreamOpenRead(cdoStreamName(0));
vlistID1 = pstreamInqVlist(streamID1);
}
else
if ( cdoVerbose ) cdoPrint("Found %i timeSteps", nts);
pstreamClose(streamID1);
streamID1 = pstreamOpenRead(cdoStreamName(0));
vlistID1 = pstreamInqVlist(streamID1);
int taxisID1 = vlistInqTaxis(vlistID1);
/* reset the requested number of eigen-function to the maximum if neccessary */
......@@ -144,20 +130,29 @@ void *EOF3d(void * argument)
if ( cdoVerbose ) cdoPrint("counted %i timesteps",n);
int nvars = vlistNvars(vlistID1);
int nrecs;
int gridID1 = vlistInqVarGrid(vlistID1, 0);
size_t gridsizemax = vlistGridsizeMax(vlistID1);
/* allocation of temporary fields and output structures */
double *in = (double *) Malloc(gridsize*sizeof(double));
double *in = (double *) Malloc(gridsizemax*sizeof(double));
int **datacounts = (int **) Malloc(nvars*sizeof(int*));
double ***datafields = (double ***) Malloc(nvars*sizeof(double **));
double ***eigenvectors = (double ***) Malloc(nvars*sizeof(double **));
double ***eigenvalues = (double ***) Malloc(nvars*sizeof(double **));
size_t maxlevs = 0;
for ( varID = 0; varID < nvars; ++varID )
{
gridsize = vlistGridsizeMax(vlistID1);
size_t gridsize = vlistGridsizeMax(vlistID1);
nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
temp_size = ((size_t)gridsize) * nlevs;
temp_size = gridsize * nlevs;
missval = vlistInqVarMissval(vlistID1, varID);
if ( nlevs > maxlevs ) maxlevs = nlevs;
datacounts[varID] = (int*) Malloc(nlevs*sizeof(int));
datafields[varID] = (double **) Malloc(nts*sizeof(double *));
......@@ -187,9 +182,28 @@ void *EOF3d(void * argument)
}
if ( cdoVerbose)
cdoPrint("allocated eigenvalue/eigenvector with nts=%i, n=%i, gridsize=%i for processing in %s",
nts,n,gridsize,"time_space");
cdoPrint("Allocated eigenvalue/eigenvector with nts=%i, n=%i, gridsize=%zu for processing in %s",
nts, n, gridsizemax, "time_space");
double *weight = (double *) Malloc(maxlevs*gridsizemax*sizeof(double));
for ( size_t i = 0; i < maxlevs*gridsizemax; ++i ) weight[i] = 1.;
if ( weight_mode == WEIGHT_ON )
{
int wstatus = gridWeights(gridID1, weight);
if ( wstatus != 0 )
{
weight_mode = WEIGHT_OFF;
cdoWarning("Using constant grid cell area weights!");
}
else
{
for ( size_t k = 1; k < maxlevs; ++k )
for ( size_t i = 0; i < gridsizemax; ++i )
weight[k*gridsizemax+i] = weight[i];
}
}
int tsID = 0;
/* read the data and create covariance matrices for each var & level */
......@@ -202,13 +216,13 @@ void *EOF3d(void * argument)
{
pstreamInqRecord(streamID1, &varID, &levelID);
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
size_t gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
missval = vlistInqVarMissval(vlistID1, varID);
pstreamReadRecord(streamID1, in, &nmiss);
offset = gridsize * levelID;
for ( int i = 0; i < gridsize; ++i )
for ( size_t i = 0; i < gridsize; ++i )
{
if ( ! DBL_IS_EQUAL(in[i], missval ) )
{
......@@ -233,13 +247,13 @@ void *EOF3d(void * argument)
if ( cdoVerbose )
cdoPrint("Read data for %i variables",nvars);
int *pack = (int*) Malloc(temp_size*sizeof(int)); //TODO
size_t *pack = (size_t*) Malloc(temp_size*sizeof(size_t)); //TODO
for ( varID = 0; varID < nvars; varID++ )
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
size_t gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
temp_size = ((size_t)gridsize) * nlevs;
temp_size = gridsize * nlevs;
if ( cdoVerbose )
{
......@@ -301,7 +315,7 @@ void *EOF3d(void * argument)
double *df2p = datafields[varID][j2];
double sum = 0;
for ( size_t i = 0; i < npack; i++ )
sum += weight[pack[i]%gridsize]*df1p[pack[i]]*df2p[pack[i]];
sum += weight[pack[i]%gridsizemax]*df1p[pack[i]]*df2p[pack[i]];
cov[j2][j1] = cov[j1][j2] = sum / sum_w / nts;
}
}
......@@ -324,7 +338,7 @@ void *EOF3d(void * argument)
/* NOW: cov contains the eigenvectors, eigv the eigenvalues */
if ( cdoVerbose )
cdoPrint("Processed SVD decomposition for var %i from %i x %i matrix",varID,n,n);
cdoPrint("Processed SVD decomposition for var %i from %zu x %zu matrix",varID,n,n);
for( int eofID = 0; eofID < n; eofID++ )
eigenvalues[varID][eofID][0] = eigv[eofID];
......@@ -351,10 +365,10 @@ void *EOF3d(void * argument)
double sum = 0;
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(eigenvec,weight,pack,npack,gridsize) reduction(+:sum)
#pragma omp parallel for default(none) shared(eigenvec,weight,pack,npack,gridsizemax) reduction(+:sum)
#endif
for ( size_t i = 0; i < npack; i++ )
sum += weight[pack[i]%gridsize] *
sum += weight[pack[i]%gridsizemax] *
eigenvec[pack[i]] * eigenvec[pack[i]];
if ( sum > 0 )
......@@ -371,7 +385,7 @@ void *EOF3d(void * argument)
#if defined(_OPENMP)
#pragma omp parallel for default(none) shared(eigenvec,pack,missval,npack)
#endif
for( size_t i = 0; i < npack; i++ )
for ( size_t i = 0; i < npack; i++ )
eigenvec[pack[i]] = missval;
}
} /* for ( eofID = 0; eofID < n_eig; eofID++ ) */
......@@ -386,7 +400,11 @@ void *EOF3d(void * argument)
/* eigenvalues */
int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
int vlistID2 = vlistDuplicate(vlistID1);
int taxisID2 = taxisDuplicate(taxisID1);
taxisDefRdate(taxisID2, 0);
taxisDefRtime(taxisID2, 0);
vlistDefTaxis(vlistID2, taxisID2);
int gridID2 = gridCreate(GRID_LONLAT, 1);
gridDefXsize(gridID2, 1);
......@@ -395,24 +413,21 @@ void *EOF3d(void * argument)
gridDefXvals(gridID2, &xvals);
gridDefYvals(gridID2, &yvals);
ngrids = vlistNgrids(vlistID2);
for ( int i = 0; i < ngrids; i++ )
vlistChangeGridIndex(vlistID2, i, gridID2);
int zaxisID2 = zaxisCreate(ZAXIS_GENERIC, 1);
double zvals = 0;
zaxisDefLevels(zaxisID2, &zvals);
zaxisDefName(zaxisID2, "zaxis_Reduced");
zaxisDefLongname(zaxisID2, "Reduced zaxis from EOF3D - only one eigen value per 3D eigen vector");
int vlistID2 = vlistCreate();
taxisDefRdate(taxisID2, 0);
taxisDefRtime(taxisID2, 0);
vlistDefTaxis(vlistID2, taxisID2);
int *varID2 = (int*) Malloc(nvars*sizeof(int));
for ( varID=0; varID<nvars; varID++ )
varID2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, TSTEP_INSTANT);
ngrids = vlistNgrids(vlistID2);
for ( int i = 0; i < ngrids; i++ )
vlistChangeGridIndex(vlistID2, i, gridID2);
int nzaxis = vlistNzaxis(vlistID2);
for ( int i = 0; i < nzaxis; i++ )
vlistChangeZaxisIndex(vlistID2, i, zaxisID2);
/* eigenvectors */
int streamID3 = pstreamOpenWrite(cdoStreamName(2), cdoFiletype());
int vlistID3 = vlistDuplicate(vlistID1);
......@@ -446,13 +461,13 @@ void *EOF3d(void * argument)
for ( varID = 0; varID < nvars; varID++ )
{
nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for ( levelID = 0; levelID < nlevs; levelID++ )
for ( levelID = 0; levelID < (int)nlevs; levelID++ )
{
offset = levelID * gridsize;
offset = levelID * gridsizemax;
if ( tsID < n_eig )
{