Eof3d.cc 14.6 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-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.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
enum T_EIGEN_MODE get_eigenmode(void);
44
45
enum T_WEIGHT_MODE get_weightmode(void);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
46

47
48
49
50
51
52
// NO MISSING VALUE SUPPORT ADDED SO FAR

void *EOF3d(void * argument)
{
  enum {EOF3D_, EOF3D_TIME, EOF3D_SPATIAL};

Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
  size_t temp_size = 0, npack = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
  int varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
  int missval_warning = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
  int nmiss, ngrids, n = 0, nlevs = 0;
57
  int offset;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
  int timer_cov = 0, timer_eig = 0;
59

60
61
62
  int calendar = CALENDAR_STANDARD;
  juldate_t juldate;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
  double missval = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
  double sum_w;
65
66
67
  double **cov = NULL;                                /* TODO: covariance matrix / eigenvectors after solving */
  double *eigv;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
69
70
71
72
  if ( cdoTimer )
    {
      timer_cov  = timer_new("Timeof cov");
      timer_eig  = timer_new("Timeof eig");
    }
73
74

  cdoInitialize(argument);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76

  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79
  cdoOperatorAdd("eof3d",        EOF3D_,        0, NULL);
  cdoOperatorAdd("eof3dtime",    EOF3D_TIME,    0, NULL);
  cdoOperatorAdd("eof3dspatial", EOF3D_SPATIAL, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
  // clang-format on
81

Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
  int operatorID  = cdoOperatorID();
  int operfunc    = cdoOperatorF1(operatorID);
84
85

  operatorInputArg("Number of eigen functions to write out");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
  int n_eig = parameter2int(operatorArgv()[0]);
87

88
  enum T_EIGEN_MODE eigen_mode = get_eigenmode();
89
  enum T_WEIGHT_MODE weight_mode = get_weightmode();
90

91
92
  int streamID1  = pstreamOpenRead(cdoStreamName(0));
  int vlistID1   = pstreamInqVlist(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
95
  int gridID1    = vlistInqVarGrid(vlistID1, 0);
  long gridsize  = vlistGridsizeMax(vlistID1);
  int nvars      = vlistNvars(vlistID1);
96
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97

98
  double *weight = (double *) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
  for ( int i = 0; i < gridsize; ++i ) weight[i] = 1.;
100

101
102
103
  if ( weight_mode == WEIGHT_ON )
    {
      int wstatus = gridWeights(gridID1, weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
      if ( wstatus != 0 )
105
106
107
108
109
	{
	  weight_mode = WEIGHT_OFF;
	  cdoWarning("Using constant grid cell area weights!");
	}
    }
110

111
  /*  eigenvalues */
112
113
114
115
116

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

  /* COUNT NUMBER OF TIMESTEPS if EOF3D_ or EOF3D_TIME */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
  int nts = vlistNtsteps(vlistID1);
118
  if ( nts == -1 )
119
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
      nts = 0;
121
      while ( pstreamInqTimestep(streamID1, nts) ) nts++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122

123
      if ( cdoVerbose ) cdoPrint("Counted %i timeSteps", nts);
124
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
126
  else
    if ( cdoVerbose ) cdoPrint("Found %i timeSteps", nts);
127

128
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129

130
131
  streamID1 = pstreamOpenRead(cdoStreamName(0));
  vlistID1  = pstreamInqVlist(streamID1);
132
  int taxisID1  = vlistInqTaxis(vlistID1);
133
134
135
136
137
138
139
140
141
142

  /* 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;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
  n = nts;
144

Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
  if ( cdoVerbose )  cdoPrint("counted %i timesteps",n);
146
147

  /* allocation of temporary fields and output structures */
148
149
150
151
152
  double *in       = (double *) Malloc(gridsize*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 **));
153
154
155
156
157

  for ( varID = 0; varID < nvars; ++varID )
    {
      gridsize            = vlistGridsizeMax(vlistID1);
      nlevs               = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
      temp_size           = ((size_t)gridsize) * nlevs;
159
160
      missval             = vlistInqVarMissval(vlistID1, varID);

161
162
      datacounts[varID]   = (int*) Malloc(nlevs*sizeof(int));
      datafields[varID]   = (double **) Malloc(nts*sizeof(double *));
163

Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
      for ( int tsID = 0; tsID < nts; tsID++ )
165
	{
166
	  datafields[varID][tsID] = (double *) Malloc(temp_size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
	  for ( size_t i = 0; i < temp_size; ++i ) datafields[varID][tsID][i] = 0;
168
	}
169
      datacounts[varID] = (int *) Malloc(temp_size*sizeof(int));	      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
      for( size_t i = 0; i < temp_size; i++) datacounts[varID][i] = 0;
171
      
172
173
      eigenvectors[varID] = (double **) Malloc(n_eig*sizeof(double *));
      eigenvalues[varID]  = (double **) Malloc(nts*sizeof(double *));
174

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
      for ( int i = 0; i < n; i++ )
176
177
178
	{
	  if ( i < n_eig )
	    {
179
	      eigenvectors[varID][i] = (double *) Malloc(temp_size*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
	      for ( size_t i2 = 0; i2 < temp_size; ++i2 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
		eigenvectors[varID][i][i2] = missval;
182
183
	    }
	  
184
	  eigenvalues[varID][i]    = (double *) Malloc(1*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
	  eigenvalues[varID][i][0] = missval;
186
187
	}
    }
188
189
190
191

  if ( cdoVerbose)
    cdoPrint("allocated eigenvalue/eigenvector with nts=%i, n=%i, gridsize=%i for processing in %s",
	     nts,n,gridsize,"time_space");
192
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
  int tsID = 0;
194
195
196
197

  /* read the data and create covariance matrices for each var & level */
  while ( TRUE )
    {
198
      nrecs = pstreamInqTimestep(streamID1, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
      if ( nrecs == 0 ) break;
200

201
      for ( int recID = 0; recID < nrecs; recID++ )
202
        {
203
          pstreamInqRecord(streamID1, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204

205
206
          gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
          missval = vlistInqVarMissval(vlistID1, varID);
208
          pstreamReadRecord(streamID1, in, &nmiss);
209
210

	  offset = gridsize * levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
	  for ( int i = 0; i < gridsize; ++i )
212
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
	      if ( ! DBL_IS_EQUAL(in[i], missval ) )
214
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
		  datafields[varID][tsID][offset + i] = in[i];
216
217
218
219
220
221
		  datacounts[varID][offset + i]++;
		}
	      else
		{
		  if ( missval_warning == 0 )
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
		      cdoWarning("Missing Value Support not checked for this Operator!");
223
224
225
		      cdoWarning("Does not work with changing locations of missing values in time.");
		      missval_warning = 1;
		    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
		  datafields[varID][tsID][i+offset] = 0;
227
228
229
230
231
232
		}
	    }
        }
      tsID++;
    }

233
234
235
  if ( cdoVerbose ) 
    cdoPrint("Read data for %i variables",nvars);
  
236
  int *pack = (int*) Malloc(temp_size*sizeof(int)); //TODO
237
238
239

  for ( varID = 0; varID < nvars; varID++ )
    {
240
      gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
      nlevs    = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
      temp_size = ((size_t)gridsize) * nlevs;
243

Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
246
247
248
249
250
      if ( cdoVerbose )
        {
          char vname[64];
          vlistInqVarName(vlistID1,varID,&vname[0]);
          cdoPrint("============================================================================");
          cdoPrint("Calculating covariance matrix and SVD for var%i (%s)",varID,vname);
        }
251

Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
      npack = 0;    // TODO already set to 0
253
254
255

      if ( cdoTimer ) timer_start(timer_cov);
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
256
      for ( size_t i = 0; i < temp_size ; i++ )
257
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
	  if ( datacounts[varID][i] > 1 )
259
260
261
262
263
	    {
	      pack[npack] = i;
	      npack++;
	    }
	}
264

265
266
267
268
      sum_w = 1;
      if ( weight_mode == WEIGHT_ON )
	{
	  sum_w = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
	  for ( size_t i = 0; i < npack; i++ )  sum_w += weight[pack[i]];
270
271
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
272
273
274
275
276
277
278
      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;
        }
279

280
	  
281
      cov = (double **) Malloc(nts*sizeof(double*));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282
      for ( int j1 = 0; j1 < nts; j1++)
283
284
	cov[j1] = (double *) Malloc(nts*sizeof(double));
      eigv = (double *) Malloc(n*sizeof(double));
285

Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
287
288
289
290
      if ( cdoVerbose )
        {
          cdoPrint("varID %i allocated eigv and cov with nts=%i and n=%i", varID, nts, n);
          cdoPrint("   npack=%zu, nts=%i temp_size=%zu", npack, nts, temp_size);
        }
291

292

Uwe Schulzweida's avatar
Uwe Schulzweida committed
293
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
#pragma omp parallel for default(shared) schedule(static,2000)
295
#endif 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
297
298
299
300
301
302
303
304
305
306
307
      for ( int j1 = 0; j1 < nts; j1++ )
        {
          double *df1p = datafields[varID][j1];
          for ( int j2 = j1; j2 < nts; j2++ )
            {
              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]];
              cov[j2][j1] = cov[j1][j2] = sum / sum_w / nts;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
309
      
      if ( cdoVerbose ) cdoPrint("calculated cov-matrix");
310

311
      /* SOLVE THE EIGEN PROBLEM */
312
      if ( cdoTimer ) timer_stop(timer_cov);
313
314


315
316
      if ( cdoTimer ) timer_start(timer_eig);

317
      if ( cdoVerbose ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
318
	cdoPrint("Processed correlation matrix for var %2i | npack: %zu", varID, n);
319
320

      if ( eigen_mode == JACOBI ) 
321
	parallel_eigen_solution_of_symmetric_matrix(cov, eigv, n, __func__);
322
      else 
323
	eigen_solution_of_symmetric_matrix(cov, eigv, n, __func__);
324
325
      /* NOW: cov contains the eigenvectors, eigv the eigenvalues */

326
327
328
      if ( cdoVerbose ) 
	cdoPrint("Processed SVD decomposition for var %i from %i x %i matrix",varID,n,n);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
      for( int eofID = 0; eofID < n; eofID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
	eigenvalues[varID][eofID][0] = eigv[eofID];
331
      
332
333
      if ( cdoTimer ) timer_stop(timer_eig);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
      for ( int eofID = 0; eofID < n_eig; eofID++ )
335
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
	  double *eigenvec = eigenvectors[varID][eofID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337

Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
#pragma omp parallel for default(none) shared(varID,nts,eofID,npack,pack,cov,datafields,eigenvec)
340
#endif 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
341
	  for ( size_t i = 0; i < npack; i++ )
342
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
344
	      double sum = 0;
	      for ( int j = 0; j < nts; j++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
346
		sum += datafields[varID][j][pack[i]] * cov[eofID][j];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
	      eigenvec[pack[i]] = sum;
348
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349

350
	  // NORMALIZING
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
	  double sum = 0;
352

Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
#if defined(_OPENMP)
354
#pragma omp parallel for default(none)  shared(eigenvec,weight,pack,npack,gridsize) reduction(+:sum)
355
#endif 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
	  for ( size_t i = 0; i < npack; i++ )
357
358
	    sum +=  weight[pack[i]%gridsize] *
	            eigenvec[pack[i]] * eigenvec[pack[i]];
359

360
361
362
	  if ( sum > 0 )
	    {
	      sum = sqrt(sum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363
#if defined(_OPENMP)
364
#pragma omp parallel for default(none) shared(sum,npack,eigenvec,pack)
365
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
	      for ( size_t i = 0; i < npack; i++ )
367
368
		eigenvec[pack[i]] /= sum;
	    }
369
	  else
370
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
#if defined(_OPENMP)
372
#pragma omp parallel for default(none) shared(eigenvec,pack,missval,npack)
373
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
	      for( size_t i = 0; i < npack; i++ )
375
376
		eigenvec[pack[i]] = missval;
	    }
377
378
	}     /* for ( eofID = 0; eofID < n_eig; eofID++ )     */

379
      if ( eigv ) Free(eigv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
      for ( int i = 0; i < n; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
381
	if ( cov[i] ) Free(cov[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
382
    }         /* for ( varID = 0; varID < nvars; varID++ )    */
383
384
385

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
  /*  eigenvalues */
387
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
388

Uwe Schulzweida's avatar
Uwe Schulzweida committed
389
  int taxisID2 = taxisDuplicate(taxisID1);
390

Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
  int gridID2 = gridCreate(GRID_LONLAT, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
393
  gridDefXsize(gridID2, 1);
  gridDefYsize(gridID2, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
395
396
397
398
399
400
401
402
  double xvals = 0, yvals = 0;
  gridDefXvals(gridID2, &xvals);
  gridDefYvals(gridID2, &yvals);

  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");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403
404
405
406
407
408

  int vlistID2 = vlistCreate();
  taxisDefRdate(taxisID2, 0);
  taxisDefRtime(taxisID2, 0);
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
  int *varID2 = (int*) Malloc(nvars*sizeof(int));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
411
  for ( varID=0; varID<nvars; varID++ )
    varID2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, TSTEP_INSTANT);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
  ngrids = vlistNgrids(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413
  for ( int i = 0; i < ngrids; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
415
    vlistChangeGridIndex(vlistID2, i, gridID2);

416
  int streamID3 = pstreamOpenWrite(cdoStreamName(2), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417

Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
419
  int vlistID3 = vlistDuplicate(vlistID1);
  int taxisID3 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
421
422
  taxisDefRdate(taxisID3, 0);
  taxisDefRtime(taxisID3, 0);
  vlistDefTaxis(vlistID3, taxisID3);
423

424
425
  pstreamDefVlist(streamID2, vlistID2);
  pstreamDefVlist(streamID3, vlistID3);
426

Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
428
  int vdate = 10101;
  int vtime = 0;
429
  juldate = juldate_encode(calendar, vdate, vtime);
430
431
  for ( tsID = 0; tsID < n; tsID++ )
    {
432
433
434
435
436
      juldate = juldate_add_seconds(60, juldate);
      juldate_decode(calendar, juldate, &vdate, &vtime);

      taxisDefVdate(taxisID2, vdate);
      taxisDefVtime(taxisID2, vtime);
437
      pstreamDefTimestep(streamID2, tsID);
438
439
440

      if ( tsID < n_eig )
        {
441
442
          taxisDefVdate(taxisID3, vdate);
          taxisDefVtime(taxisID3, vtime);
443
          pstreamDefTimestep(streamID3, tsID);
444
445
446
447
448
449
450
451
452
453
454
        }

      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;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
                  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
                    if ( DBL_IS_EQUAL(eigenvectors[varID][tsID][offset + i], missval) ) nmiss++;
457

458
459
                  pstreamDefRecord(streamID3, varID, levelID);
                  pstreamWriteRecord(streamID3, &eigenvectors[varID][tsID][offset], nmiss);
460
461
                }
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
462
	  if ( DBL_IS_EQUAL(eigenvalues[varID][tsID][0], missval) ) nmiss = 1;
463
	  else nmiss = 0;
464
465
	  pstreamDefRecord(streamID2, varID, 0);
	  pstreamWriteRecord(streamID2, eigenvalues[varID][tsID],nmiss);
466
467
468
469
470
        } // for ( varID = 0; ... )
    } // for ( tsID = 0; ... )

  for ( varID = 0; varID < nvars; varID++)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
      for ( int i = 0; i < nts; i++)
472
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
	  Free(datafields[varID][i]);
474
	  if ( i < n_eig )
475
476
	    Free(eigenvectors[varID][i]);
	  Free(eigenvalues[varID][i]);
477
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478

479
480
481
482
      Free(datafields[varID]);
      Free(datacounts[varID]);
      Free(eigenvectors[varID]);
      Free(eigenvalues[varID]);
483
484
    }

485
486
487
488
489
  Free(datafields);
  Free(datacounts);
  Free(eigenvectors);
  Free(eigenvalues);
  Free(in);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
490
  Free(varID2);
491

492
493
  Free(pack);
  Free(weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
494

495
496
497
  pstreamClose(streamID3);
  pstreamClose(streamID2);
  pstreamClose(streamID1);
498
499
500
  
  cdoFinish();
 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
  return 0;
502
}