Eof3d.c 17.1 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2014 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

/*
   This module contains the following operators:

     EOF3d        eof3d             3D-EOF in spatial or time space
     EOF3d        eof3dspatial      3D-EOF in spatial space
     EOF3d        eof3dtime         3D-EOF in time space
*/
/*
 * TODO: 
 * Role of the weights for eofs. Should not be mixed up with division with
 * number of contributing values during summation.
 */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
#include <omp.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
#endif
34
35
36
37
38

#include "cdi.h"
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
#include "grid.h"
40
#include "statistic.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
43
44
45


enum T_EIGEN_MODE {JACOBI, DANIELSON_LANCZOS};

#define WEIGHTS 1
46
47
48
49
50
51
52
53
54
55

// NO MISSING VALUE SUPPORT ADDED SO FAR

void *EOF3d(void * argument)
{
  char *envstr;

  enum {EOF3D_, EOF3D_TIME, EOF3D_SPATIAL};

  int **datacounts;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
  int gridsize, temp_size = 0;
57
58
59
60
61
62
63
64
65
66
  int gridID1, gridID2, gridID3;
  int i, i2, j, j1, j2, eofID, varID, recID, levelID, tsID;
  int missval_warning=0;
  int nmiss,ngrids,n_eig,nrecs,nvars,n=0,nlevs=0,npack=0,nts=0;
  int offset;
  int operatorID, operfunc;
  int *pack;
  int reached_eof;
  int streamID1, streamID2, streamID3;
  int taxisID1, taxisID2, taxisID3;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
  int timer_init = 0, timer_alloc = 0, timer_read = 0, timer_cov = 0;
  int timer_eig = 0, timer_post = 0, timer_write = 0, timer_finish = 0;
69
70
  int *varID2;
  int vdate=0, vtime=0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
  int vlistID1, vlistID2 = -1, vlistID3 = -1;
72
73
  int zaxisID2;

74
75
76
  int calendar = CALENDAR_STANDARD;
  juldate_t juldate;

77
78
79
80
81
82
  double missval=0;
  double sum_w, sum;
  double **cov = NULL;                                /* TODO: covariance matrix / eigenvectors after solving */
  double *eigv;
  double *weight;
  double *xvals, *yvals, *zvals;
83
  double *df1p, *df2p;
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

  field_t **datafields;
  field_t **eigenvectors, **eigenvalues;
  field_t in;

  enum T_EIGEN_MODE eigen_mode = JACOBI;


  if ( cdoTimer ) {
    timer_init = timer_new("Timeof init");
    timer_alloc= timer_new("Timeof alloc");
    timer_read = timer_new("Timeof read");
    timer_cov  = timer_new("Timeof cov");
    timer_eig  = timer_new("Timeof eig");
    timer_post = timer_new("Timeof post");
    timer_write= timer_new("Timeof write");
    timer_finish=timer_new("Timeof finish");

    timer_start(timer_init);
  }

  cdoInitialize(argument);
  cdoOperatorAdd("eof3d",       EOF3D_,       0, NULL);
  cdoOperatorAdd("eof3dtime",   EOF3D_TIME,   0, NULL);
  cdoOperatorAdd("eof3dspatial",EOF3D_SPATIAL,0, NULL);

  operatorID  = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
  operfunc    = cdoOperatorF1(operatorID);
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

  operatorInputArg("Number of eigen functions to write out");
  n_eig       = atoi(operatorArgv()[0]);

  envstr = getenv("CDO_SVD_MODE");

  if ( envstr &&! strncmp(envstr,"danielson_lanczos",17) )
    eigen_mode = DANIELSON_LANCZOS;
  else if ( envstr && ! strncmp(envstr,"jacobi",6) )
    eigen_mode = JACOBI;
  else if ( envstr ) {
    cdoWarning("Unknown environmental setting %s for CDO_SVD_MODE. Available options are",envstr);
    cdoWarning("  - 'jacobi' for a one-sided parallelized jacobi algorithm");
    cdoWarning("  - 'danielson_lanzcos' for the D/L algorithm");
  }

  if ( cdoVerbose ) 
    cdoPrint("Set eigen_mode to %s\n",eigen_mode == JACOBI? "jacobi" : "danielson_lanczos");

Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
#if defined(_OPENMP)
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
  if ( omp_get_max_threads() > 1 && eigen_mode == DANIELSON_LANCZOS )  {
    cdoWarning("Requested parallel computation with %i Threads ",omp_get_max_threads());
    cdoWarning("  but environmental setting CDO_SVD_MODE causes sequential ");
    cdoWarning("  Singular value decomposition");
  }
#endif 

  streamID1   = streamOpenRead(cdoStreamName(0));
  vlistID1    = streamInqVlist(streamID1);
  taxisID1    = vlistInqTaxis(vlistID1);
  gridID1     = vlistInqVarGrid(vlistID1, 0);
  gridsize    = vlistGridsizeMax(vlistID1);
  nvars       = vlistNvars(vlistID1);
  nrecs       = vlistNrecs(vlistID1);
  taxisID1    = vlistInqTaxis(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
  weight      = (double*) malloc(gridsize*sizeof(double));
148
149
150
151
152
153
154
155
156
  if ( WEIGHTS )
      gridWeights(gridID1, &weight[0]);
  else
    for(i=0;i<gridsize;i++)
      weight[i]=1;


  /*  eigenvalues */
  streamID2   = streamOpenWrite(cdoStreamName(1), cdoFiletype());
157

158
159
160
161
162
  taxisID2    = taxisDuplicate(taxisID1);

  gridID2     = gridCreate(GRID_LONLAT, 1);
  gridDefXsize(gridID2, 1);
  gridDefYsize(gridID2, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
164
165
  xvals       = (double*) malloc(1*sizeof(double));
  yvals       = (double*) malloc(1*sizeof(double));
  zvals       = (double*) malloc(1*sizeof(double));
166
167
168
169
170
171
172
173
174
175
176
177
  xvals[0]    = 0;
  yvals[0]    = 0;
  zvals[0]    = 0;
  gridDefXvals(gridID2, xvals);
  gridDefYvals(gridID2, yvals);

  zaxisID2 = zaxisCreate(ZAXIS_GENERIC,1);
  zaxisDefLevels(zaxisID2,zvals);
  zaxisDefName(zaxisID2,"zaxis_Reduced");
  zaxisDefLongname(zaxisID2,"Reduced zaxis from EOF3D - only one eigen value per 3D eigen vector");

  vlistID2 = vlistCreate();
178
179
  taxisDefRdate(taxisID2, 0);
  taxisDefRtime(taxisID2, 0);
180
181
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
  varID2 = (int*) malloc(nvars*sizeof(int));
183
  for ( varID=0; varID<nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
    varID2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, TSTEP_INSTANT);
185
186
187
188
189
190
  ngrids      = vlistNgrids(vlistID2);
  for ( i = 0; i < ngrids; i++ )
    vlistChangeGridIndex(vlistID2, i, gridID2);

  /*  eigenvectors */
  streamID3   = streamOpenWrite(cdoStreamName(2), cdoFiletype());
191

192
193
194
  vlistID3    = vlistDuplicate(vlistID1);
  taxisID3    = taxisDuplicate(taxisID1);
  gridID3     = gridDuplicate(gridID1);
195
196
  taxisDefRdate(taxisID3, 0);
  taxisDefRtime(taxisID3, 0);
197
198
199
  vlistDefTaxis(vlistID3, taxisID3);
  

200
201
202
  if ( cdoVerbose )
    cdoPrint("Initialized streams");

203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
  /*  eigenvalues */

  reached_eof = 0;
  tsID        = 0;

  if ( operfunc == EOF3D_SPATIAL )
    cdoAbort("Operator not Implemented - use eof3d or eof3dtime instead");


  /* COUNT NUMBER OF TIMESTEPS if EOF3D_ or EOF3D_TIME */
  while ( TRUE )
    {
      if ( reached_eof ) continue;
      nrecs = streamInqTimestep(streamID1, tsID);
      if ( nrecs == 0 ) {
	reached_eof = 1;
	break;
      }
      tsID++;
    }

  nts         = tsID;
  reached_eof = 0;
  streamID1   = streamOpenRead(cdoStreamName(0));

  /* reset the requested number of eigen-function to the maximum if neccessary */
  if ( n_eig > nts )
    {
      cdoWarning("Solving in time-space:");
      cdoWarning("Number of eigen-functions to write out is bigger than number of time-steps.");
      cdoWarning("Setting n_eig to %i.", nts);
      n_eig = nts;
    }
  n = nts;

238
239
240
  if ( cdoVerbose ) 
    cdoPrint("counted %i timesteps",n);

241
242
243
244
  if ( cdoTimer ) timer_stop(timer_init);
  if ( cdoTimer ) timer_start(timer_alloc);

  /* allocation of temporary fields and output structures */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
246
  in.ptr       = (double*) malloc(gridsize*sizeof(double));
  datafields   = (field_t**) malloc(nvars*sizeof(field_t*));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
  datacounts   = (int**) malloc(nvars*sizeof(int*));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
  eigenvectors = (field_t**) malloc(nvars*sizeof(field_t*));
  eigenvalues  = (field_t**) malloc(nvars*sizeof(field_t*));
250
251
252
253
254
255
256
257
258

  for ( varID = 0; varID < nvars; ++varID )
    {
      gridID1             = vlistInqVarGrid(vlistID1, varID);
      gridsize            = vlistGridsizeMax(vlistID1);
      nlevs               = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
      temp_size           = gridsize * nlevs;
      missval             = vlistInqVarMissval(vlistID1, varID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
260
261
      datacounts[varID]   = (int*) malloc(nlevs*sizeof(int));
      eigenvectors[varID] = (field_t*) malloc(nlevs*sizeof(field_t));
      datafields[varID]   = (field_t*) malloc(nts*sizeof(field_t));
262
263
264
265
266
267

      for ( tsID = 0; tsID < nts; tsID++ )
	{
	  datafields[varID][tsID].grid    = gridID1;
	  datafields[varID][tsID].nmiss   = 0;
	  datafields[varID][tsID].missval = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
268
	  datafields[varID][tsID].ptr     = (double*) malloc(temp_size*sizeof(double));
269
270
271
	  for ( i = 0; i < temp_size; ++i )
	    datafields[varID][tsID].ptr[i] = 0;
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
272
      datacounts[varID] = (int*) malloc(temp_size*sizeof(int));	      
273
274
275
      for(i=0;i<temp_size;i++)
	datacounts[varID][i] = 0;
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276
277
      eigenvectors[varID] = (field_t*) malloc(n_eig*sizeof(field_t));
      eigenvalues[varID]  = (field_t*) malloc(nts*sizeof(field_t));
278

279
280
281
282
283
284
285
      for ( i = 0; i < n; i++ )
	{
	  if ( i < n_eig )
	    {
	      eigenvectors[varID][i].grid    = gridID2;
	      eigenvectors[varID][i].nmiss   = 0;
	      eigenvectors[varID][i].missval = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
	      eigenvectors[varID][i].ptr     = (double*) malloc(temp_size*sizeof(double));
287
288
289
290
291
292
293
	      for ( i2 = 0; i2 < temp_size; ++i2 )
		eigenvectors[varID][i].ptr[i2] = missval;
	    }
	  
	  eigenvalues[varID][i].grid    = gridID3;
	  eigenvalues[varID][i].nmiss   = 0;
	  eigenvalues[varID][i].missval = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
	  eigenvalues[varID][i].ptr     = (double*) malloc(1*sizeof(double));
295
296
297
	  eigenvalues[varID][i].ptr[0]  = missval;
	}
    }
298
299
300
301

  if ( cdoVerbose)
    cdoPrint("allocated eigenvalue/eigenvector with nts=%i, n=%i, gridsize=%i for processing in %s",
	     nts,n,gridsize,"time_space");
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
  
  if ( cdoTimer ) timer_stop(timer_alloc);

  if ( cdoTimer ) timer_start(timer_read);
  tsID = 0;

  /* read the data and create covariance matrices for each var & level */
  while ( TRUE )
    {
      if ( reached_eof ) continue;
      nrecs = streamInqTimestep(streamID1, tsID);
      if ( nrecs == 0 )
        {
          reached_eof = 1;
          break;
        }

      vdate = taxisInqVdate(taxisID1);
      vtime = taxisInqVtime(taxisID1);

      for ( recID = 0; recID < nrecs; recID++ )
        {
          streamInqRecord(streamID1, &varID, &levelID);
          gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));

          missval  = in.missval = vlistInqVarMissval(vlistID1, varID);
          streamReadRecord(streamID1, in.ptr, &in.nmiss);

	  offset = gridsize * levelID;
	  for ( i=0; i<gridsize; ++i )
	    {
	      if ( ! DBL_IS_EQUAL(in.ptr[i], missval ) )
		{
		  datafields[varID][tsID].ptr[offset + i] = in.ptr[i];
		  datacounts[varID][offset + i]++;
		}
	      else
		{
		  if ( missval_warning == 0 )
		    {
		      cdoWarning("Missing Value Support not Checked for this Operator!");
		      cdoWarning("Does not work with changing locations of missing values in time.");
		      missval_warning = 1;
		    }
		  datafields[varID][tsID].ptr[i+offset] = 0;
		}
	    }
        }
      tsID++;
    }

353
354
355
  if ( cdoVerbose ) 
    cdoPrint("Read data for %i variables",nvars);
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
  pack = (int*) malloc(temp_size*sizeof(int)); //TODO
357
358
359
360
361

  if ( cdoTimer ) timer_stop(timer_read);

  for ( varID = 0; varID < nvars; varID++ )
    {
362
363
364
365
366
367
368
369
370
371
372
      gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
      nlevs               = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
      temp_size = gridsize * nlevs;

      if ( cdoVerbose )  {
	char vname[64];
	vlistInqVarName(vlistID1,varID,&vname[0]);
	cdoPrint("============================================================================");
	cdoPrint("Calculating covariance matrix and SVD for var%i (%s)",varID,vname);
      }

373
374
375
376
377
      npack        = 0;    // TODO already set to 0
      sum_w        = 0;

      if ( cdoTimer ) timer_start(timer_cov);
      
378
379
      // sum_w = 0;
      sum_w = 1;
380
381
      for ( i = 0; i < temp_size ; i++ )
	{
382
	  if ( datacounts[varID][i] > 1)
383
384
385
	    {
	      pack[npack] = i;
	      npack++;
386
	      //  sum_w += weight[i%gridsize];
387
388
	    }
	}
389
390
391
392
393
394
395
396

      if ( npack < 1 ) {
	char vname[64];
	vlistInqVarName(vlistID1,varID,&vname[0]);
	cdoWarning("Refusing to calculate EOF from a single time step for var%i (%s)",varID,&vname[0]);
	continue;
      }

397
	  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
      cov = (double**) malloc(nts*sizeof(double*));
399
      for ( j1 = 0; j1 < nts; j1++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
401
	cov[j1] = (double*) malloc(nts*sizeof(double));
      eigv = (double*) malloc(n*sizeof(double));
402
403
404
405
406
407

      if ( cdoVerbose )  {
	cdoPrint("varID %i allocated eigv and cov with nts=%i and n=%i",varID,nts,n);
	cdoPrint("   npack=%i, nts=%i temp_size=%i",npack,nts,temp_size);
      }

408

Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
#if defined(_OPENMP)
410
411
#pragma omp parallel for private(j1,j2,sum,df1p,df2p) default(shared) schedule(static,2000)
#endif 
412
413
414
415
      for ( j1 = 0; j1 < nts; j1++)
	for ( j2 = j1; j2 < nts; j2++ )
	  {
	    sum = 0;
416
417
	    df1p = datafields[varID][j1].ptr;
	    df2p = datafields[varID][j2].ptr;
418
	    for ( i = 0; i < npack; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
420
	      //  sum += df1p[pack[i]]*df2p[pack[i]];
	      sum += weight[pack[i]%gridsize]*df1p[pack[i]]*df2p[pack[i]];
421
422
	    cov[j2][j1] = cov[j1][j2] = sum / sum_w / nts;
	  }
423
424
      if ( cdoVerbose ) 
	cdoPrint("calculated cov-matrix");
425

426
      /* SOLVE THE EIGEN PROBLEM */
427
      if ( cdoTimer ) timer_stop(timer_cov);
428
429


430
431
      if ( cdoTimer ) timer_start(timer_eig);

432
433
      if ( cdoVerbose ) 
	cdoPrint("Processed correlation matrix for var %2i | npack: %4i",varID,n);
434
435

      if ( eigen_mode == JACOBI ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
	parallel_eigen_solution_of_symmetric_matrix(&cov[0],&eigv[0],n,n,__func__);
437
      else 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
	eigen_solution_of_symmetric_matrix(&cov[0],&eigv[0],n,n,__func__);
439
440
      /* NOW: cov contains the eigenvectors, eigv the eigenvalues */

441
442
443
      if ( cdoVerbose ) 
	cdoPrint("Processed SVD decomposition for var %i from %i x %i matrix",varID,n,n);

444
445
      for( eofID=0; eofID<n; eofID++ )
	eigenvalues[varID][eofID].ptr[0] = eigv[eofID];
446
      
447
448
449
      if ( cdoTimer ) timer_stop(timer_eig);
      if ( cdoTimer ) timer_start(timer_post);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
450
      for ( eofID = 0; eofID < n_eig; eofID++ )
451
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452
453
	  double *eigenvec = eigenvectors[varID][eofID].ptr;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
#pragma omp parallel for private(i,j,sum) shared(datafields, eigenvec)
456
#endif 
457
458
459
460
461
	  for ( i = 0; i < npack; i++ )
	    {
	      sum = 0;
	      for ( j = 0; j < nts; j++ )
		sum += datafields[varID][j].ptr[pack[i]] * cov[eofID][j];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
	      eigenvec[pack[i]] = sum;
463
464
465
	    }
	  // NORMALIZING
	  sum = 0;
466

Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
#if defined(_OPENMP)
468
#pragma omp parallel for private(i) default(none) reduction(+:sum) \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
  shared(eigenvec,weight,pack,npack,gridsize)
470
#endif 
471
	  for ( i = 0; i < npack; i++ )
472
	    // sum +=  weight[pack[i]%gridsize] *
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
	    sum += eigenvec[pack[i]] * eigenvec[pack[i]];
474
475
476

	  if ( sum > 0 ) {
	    sum = sqrt(sum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
#if defined(_OPENMP)
478
#pragma omp parallel for private(i) default(none) \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
  shared(sum,npack,eigenvec,pack)
480
#endif
481
	    for( i = 0; i < npack; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
	      eigenvec[pack[i]] /= sum;
483
	  }
484
	  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
#if defined(_OPENMP)
486
#pragma omp parallel for private(i) default(none) \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
487
  shared(eigenvec,pack,missval,npack)
488
#endif
489
	    for( i = 0; i < npack; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490
	      eigenvec[pack[i]] = missval;
491
492
493
494
	}     /* for ( eofID = 0; eofID < n_eig; eofID++ )     */

      if ( cdoTimer ) timer_stop(timer_post);

495
496
497
498
      if ( eigv ) free(eigv);
      for ( i=0; i<n; i++ )
	if ( cov[i] ) 
	  free(cov[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
499
    }         /* for ( varID = 0; varID < nvars; varID++ )    */
500
501
502

  /* write files with eigenvalues (ID3) and eigenvectors (ID2) */

503

504
  if ( cdoTimer ) timer_start(timer_write);
505
506

  cdoPrint("Started writing");
507
508
509
  streamDefVlist(streamID2, vlistID2);
  streamDefVlist(streamID3, vlistID3);

510
511
512
  vdate = 10101;
  vtime = 0;
  juldate = juldate_encode(calendar, vdate, vtime);
513
514
  for ( tsID = 0; tsID < n; tsID++ )
    {
515
516
517
518
519
      juldate = juldate_add_seconds(60, juldate);
      juldate_decode(calendar, juldate, &vdate, &vtime);

      taxisDefVdate(taxisID2, vdate);
      taxisDefVtime(taxisID2, vtime);
520
521
522
523
      streamDefTimestep(streamID2, tsID);

      if ( tsID < n_eig )
        {
524
525
          taxisDefVdate(taxisID3, vdate);
          taxisDefVtime(taxisID3, vtime);
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
          streamDefTimestep(streamID3, tsID);
        }

      for ( varID = 0; varID < nvars; varID++ )
        {
          nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
          for ( levelID = 0; levelID < nlevs; levelID++ )
            {
	      offset = levelID * gridsize;
              if ( tsID < n_eig )
                {
                  nmiss = 0;
                  for ( i = 0; i < gridsize; i++ )
                    if ( DBL_IS_EQUAL(eigenvectors[varID][tsID].ptr[offset + i], missval) ) nmiss++;

                  streamDefRecord(streamID3, varID, levelID);
                  streamWriteRecord(streamID3, &eigenvectors[varID][tsID].ptr[offset], nmiss);
                }
	    }
	  if ( DBL_IS_EQUAL(eigenvalues[varID][tsID].ptr[i], missval) ) nmiss = 1;
	  else nmiss = 0;
	  streamDefRecord(streamID2, varID, 0);
	  streamWriteRecord(streamID2, eigenvalues[varID][tsID].ptr,nmiss);
        } // for ( varID = 0; ... )
    } // for ( tsID = 0; ... )

  if ( cdoTimer ) timer_stop(timer_write);
  
  if ( cdoTimer ) timer_start(timer_finish);
555
556

  cdoPrint("Started cleanup in eof3d");
557
558
559
  
  for ( varID = 0; varID < nvars; varID++)
    {
560
      for(i = 0; i < nts; i++)
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
	{
	  if ( i < n_eig )
	    free(eigenvectors[varID][i].ptr);
	  free(eigenvalues[varID][i].ptr);
	}
      free(eigenvectors[varID]);
      free(eigenvalues[varID]);
      free(datacounts[varID]);
    }

  free(eigenvectors);
  free(eigenvalues);
  free(datafields);
  free(datacounts);
  free(in.ptr);

  streamClose(streamID3);
  streamClose(streamID2);
  streamClose(streamID1);

  if ( cdoTimer ) timer_stop(timer_finish);
  
  cdoFinish();
 
  return (0);
}