Vertstat.cc 14 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
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>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
  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.
*/

18
19
20
/*
   This module contains the following operators:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
      Vertstat   vertrange       Vertical range
22
23
24
      Vertstat   vertmin         Vertical minimum
      Vertstat   vertmax         Vertical maximum
      Vertstat   vertsum         Vertical sum
25
      Vertstat   vertint         Vertical integral
26
27
      Vertstat   vertmean        Vertical mean
      Vertstat   vertavg         Vertical average
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
      Vertstat   vertvar         Vertical variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
      Vertstat   vertvar1        Vertical variance [Normalize by (n-1)]
30
      Vertstat   vertstd         Vertical standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
      Vertstat   vertstd1        Vertical standard deviation [Normalize by (n-1)]
32
33
34
*/


Ralf Mueller's avatar
Ralf Mueller committed
35
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
38
39
40
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
#define IS_SURFACE_LEVEL(zaxisID)  (zaxisInqType(zaxisID) == ZAXIS_SURFACE && zaxisInqSize(zaxisID) == 1)

43

Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
45
46
47
48
49
50
int getSurfaceID(int vlistID)
{
  int surfID = -1;
  int nzaxis = vlistNzaxis(vlistID);

  for ( int index = 0; index < nzaxis; ++index )
    {
51
      int zaxisID = vlistZaxis(vlistID, index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
      if ( IS_SURFACE_LEVEL(zaxisID) )
	{
	  surfID = vlistZaxis(vlistID, index);
	  break;
	}
    }

  if ( surfID == -1 ) surfID = zaxisCreate(ZAXIS_SURFACE, 1);

  return surfID;
}

static
void setSurfaceID(int vlistID, int surfID)
{
  int nzaxis = vlistNzaxis(vlistID);

  for ( int index = 0; index < nzaxis; ++index )
    {
71
      int zaxisID = vlistZaxis(vlistID, index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
73
74
75
76
      if ( zaxisID != surfID || !IS_SURFACE_LEVEL(zaxisID) )
	vlistChangeZaxisIndex(vlistID, index, surfID);
    }
}

77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

void genLayerBounds(int nlev, double *levels, double *lbounds, double *ubounds)
{
  if ( nlev == 1 )
    {
      lbounds[0] = 0.;
      ubounds[0] = 1.;
    }
  else
    {
      lbounds[0]      = levels[0];
      ubounds[nlev-1] = levels[nlev-1];
      for ( int i = 0; i < nlev-1; ++i )
        {
          double bound = 0.5*(levels[i] + levels[i+1]);
          lbounds[i+1] = bound;
          ubounds[i]   = bound;
        }
    }
}

98

99
int getLayerThickness(bool genbounds, int index, int zaxisID, int nlev, double *thickness, double *weights)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
{
101
  int status = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
  int i;
103
104
105
  double *levels  = (double *) Malloc(nlev*sizeof(double));
  double *lbounds = (double *) Malloc(nlev*sizeof(double));
  double *ubounds = (double *) Malloc(nlev*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106

107
  cdoZaxisInqLevels(zaxisID, levels);
108
  if ( genbounds )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
    {
110
      status = 2;
111
      genLayerBounds(nlev, levels, lbounds, ubounds);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
    }
113
114
115
116
117
118
  else if ( zaxisInqLbounds(zaxisID, NULL) && zaxisInqUbounds(zaxisID, NULL) )
    {
      status = 1;
      zaxisInqLbounds(zaxisID, lbounds);
      zaxisInqUbounds(zaxisID, ubounds);
    }
119
120
121
122
123
124
125
126
127
  else
    {
      for ( i = 0; i < nlev; ++i )
	{
	  lbounds[i] = 0.;
	  ubounds[i] = 1.;
	}
    }
  
128
  for ( i = 0; i < nlev; ++i ) thickness[i] = fabs(ubounds[i]-lbounds[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129

130
131
132
  double lsum = 0;
  for ( i = 0; i < nlev; ++i ) lsum += thickness[i];

133
134
  for ( i = 0; i < nlev; ++i ) weights[i] = thickness[i];
  
135
136
  for ( i = 0; i < nlev; ++i ) weights[i] /= (lsum/nlev);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
  double wsum = 0;
  for ( i = 0; i < nlev; ++i ) wsum += weights[i];
139

Uwe Schulzweida's avatar
Uwe Schulzweida committed
140
141
  if ( cdoVerbose )
    {
142
      cdoPrint("zaxisID=%d  nlev=%d  layersum=%g  weightsum=%g", index, nlev, lsum, wsum);
143
      printf("         level     bounds   thickness  weight\n");
144
      for ( i = 0; i < nlev; ++i )
145
	printf("   %3d  %6g  %6g/%-6g  %6g  %6g\n", i+1, levels[i], lbounds[i], ubounds[i], thickness[i], weights[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
147
    }

148
149
150
  Free(levels);
  Free(lbounds);
  Free(ubounds);
151
152

  return status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154

155

Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
void *Vertstat(void *argument)
{
158
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
  int gridID;
160
161
  int varID, levelID;
  int nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
163
  typedef struct {
    int zaxisID;
164
    int status;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
    int numlevel;
166
    double *thickness;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
168
169
    double *weights;
  }
  vert_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
171
172

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
174
175
176
177
178
179
180
181
182
183
                 cdoOperatorAdd("vertrange", func_range, 0, NULL);
                 cdoOperatorAdd("vertmin",   func_min,   0, NULL);
                 cdoOperatorAdd("vertmax",   func_max,   0, NULL);
                 cdoOperatorAdd("vertsum",   func_sum,   0, NULL);
  int VERTINT  = cdoOperatorAdd("vertint",   func_sum,   1, NULL);
                 cdoOperatorAdd("vertmean",  func_mean,  1, NULL);
                 cdoOperatorAdd("vertavg",   func_avg,   1, NULL);
                 cdoOperatorAdd("vertvar",   func_var,   1, NULL);
                 cdoOperatorAdd("vertvar1",  func_var1,  1, NULL);
                 cdoOperatorAdd("vertstd",   func_std,   1, NULL);
                 cdoOperatorAdd("vertstd1",  func_std1,  1, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
186
187
  int operatorID   = cdoOperatorID();
  int operfunc     = cdoOperatorF1(operatorID);
  bool needWeights = cdoOperatorF2(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188

Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
  bool lrange  = operfunc == func_range;
190
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
191
192
  bool lstd    = operfunc == func_std || operfunc == func_std1;
  bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
  int  divisor = operfunc == func_std1 || operfunc == func_var1;
194

195
  //int applyWeights = lmean;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196

197
  int streamID1 = streamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

199
  int vlistID1 = streamInqVlist(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
201

  vlistClearFlag(vlistID1);
202
  int nvars = vlistNvars(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
205
  for ( varID = 0; varID < nvars; varID++ )
    vlistDefFlag(vlistID1, varID, 0, TRUE);

206
  int vlistID2 = vlistCreate();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
  vlistCopyFlag(vlistID2, vlistID1);

209
210
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
212
  vlistDefTaxis(vlistID2, taxisID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
215
216
217
  int surfID = getSurfaceID(vlistID1);
  setSurfaceID(vlistID2, surfID);

  int nzaxis = vlistNzaxis(vlistID1);
  int nlev, zaxisID;
218
  vert_t *vert = (vert_t*) Malloc(nzaxis*sizeof(vert_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
  if ( needWeights )
    {
221
      bool genbounds = false;
222
      unsigned npar = operatorArgc();
223
224
225
226
227
      if ( npar > 0 )
	{
	  char **parnames = operatorArgv();

	  if ( cdoVerbose )
228
	    for ( unsigned i = 0; i < npar; i++ )
229
230
	      cdoPrint("key %d = %s", i+1, parnames[i]);

231
	  if ( strcmp(parnames[0], "genbounds") == 0 ) genbounds = true;
232
	  else cdoAbort("Parameter >%s< unsupported! Supported parameter are: genbounds", parnames[0]);
233
234
	}
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237
238
239
      for ( int index = 0; index < nzaxis; ++index )
	{
	  zaxisID = vlistZaxis(vlistID1, index);
	  nlev = zaxisInqSize(zaxisID);
	  vert[index].numlevel = 0;
240
	  vert[index].status   = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
	  vert[index].zaxisID  = zaxisID;
242
	  // if ( nlev > 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
	    {
	      vert[index].numlevel = nlev;
245
246
	      vert[index].thickness = (double *) Malloc(nlev*sizeof(double));
	      vert[index].weights = (double *) Malloc(nlev*sizeof(double));
247
	      vert[index].status = getLayerThickness(genbounds, index, zaxisID, nlev, vert[index].thickness, vert[index].weights); 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
	    }
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251

252
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253
254
255

  streamDefVlist(streamID2, vlistID2);

256
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
257

258
  field_type field;
259
  field_init(&field);
260
  field.ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261

262
263
264
  field_type *vars1 = (field_type*) Malloc(nvars*sizeof(field_type));
  field_type *samp1 = (field_type*) Malloc(nvars*sizeof(field_type));
  field_type *vars2 = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
  if ( lvarstd || lrange )
266
    vars2 = (field_type*) Malloc(nvars*sizeof(field_type));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267
268
269
270
271

  for ( varID = 0; varID < nvars; varID++ )
    {
      gridID   = vlistInqVarGrid(vlistID1, varID);
      gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
272
      zaxisID  = vlistInqVarZaxis(vlistID1, varID);
273
      double missval = vlistInqVarMissval(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

275
276
      field_init(&vars1[varID]);
      field_init(&samp1[varID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277
      vars1[varID].grid    = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
278
      vars1[varID].zaxis   = zaxisID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
280
281
      vars1[varID].nsamp   = 0;
      vars1[varID].nmiss   = 0;
      vars1[varID].missval = missval;
282
      vars1[varID].ptr     = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
284
285
286
      samp1[varID].grid    = gridID;
      samp1[varID].nmiss   = 0;
      samp1[varID].missval = missval;
      samp1[varID].ptr     = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
287
      if ( lvarstd || lrange )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
	{
289
	  field_init(&vars2[varID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
291
292
	  vars2[varID].grid    = gridID;
	  vars2[varID].nmiss   = 0;
	  vars2[varID].missval = missval;
293
	  vars2[varID].ptr     = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
295
296
	}
    }

297
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
299
300
301
302
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);
      streamDefTimestep(streamID2, tsID);

303
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
305
	{
	  streamInqRecord(streamID1, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306

Uwe Schulzweida's avatar
Uwe Schulzweida committed
307
          vars1[varID].nsamp++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
          if ( lrange ) vars2[varID].nsamp++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
	  gridsize = gridInqSize(vars1[varID].grid);
310
311
	  zaxisID  = vars1[varID].zaxis;
	  nlev = zaxisInqSize(zaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312

313
314
315
	  double layer_weight = 1.0;
	  double layer_thickness = 1.0;
	  if ( needWeights )
316
317
318
	    {
	      for ( int index = 0; index < nzaxis; ++index )
		if ( vert[index].zaxisID == zaxisID )
319
		  {		    
320
		    if ( vert[index].status == 0 && tsID == 0 && levelID == 0 && nlev > 1 )
321
322
323
		      {
			char varname[CDI_MAX_NAME];
			vlistInqVarName(vlistID1, varID, varname);
324
			cdoWarning("Layer bounds not available, using constant vertical weights for variable %s!", varname);
325
		      }
326
327
		    else
		      {
328
329
			layer_weight    = vert[index].weights[levelID];
			layer_thickness = vert[index].thickness[levelID];
330
		      }
331
332
333
334

		    break;
		  }
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335

Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
337
338
	  if ( levelID == 0 )
	    {
	      streamReadRecord(streamID1, vars1[varID].ptr, &nmiss);
339
	      vars1[varID].nmiss = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
341
342
343
344
345
              if ( lrange )
                {
                  vars2[varID].nmiss = (size_t)nmiss;
                  for ( int i = 0; i < gridsize; i++ )
                    vars2[varID].ptr[i] = vars1[varID].ptr[i];
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
346

347
	      if ( operatorID == VERTINT && IS_NOT_EQUAL(layer_thickness, 1.0) ) farcmul(&vars1[varID], layer_thickness);
348
	      if ( lmean && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&vars1[varID], layer_weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349

350
	      if ( lvarstd )
351
352
353
354
355
356
357
358
359
360
361
                {
                  if ( IS_NOT_EQUAL(layer_weight, 1.0) )
                    {
                      farmoqw(&vars2[varID], vars1[varID], layer_weight);
                      farcmul(&vars1[varID], layer_weight);
                    }
                  else
                    {
                      farmoq(&vars2[varID], vars1[varID]);
                    }
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
362

363
	      if ( nmiss > 0 || samp1[varID].ptr || needWeights )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
364
365
		{
		  if ( samp1[varID].ptr == NULL )
366
		    samp1[varID].ptr = (double *) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
367

368
		  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
		    if ( DBL_IS_EQUAL(vars1[varID].ptr[i], vars1[varID].missval) )
370
		      samp1[varID].ptr[i] = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
		    else
372
		      samp1[varID].ptr[i] = layer_weight;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
374
375
376
		}
	    }
	  else
	    {
377
378
	      streamReadRecord(streamID1, field.ptr, &nmiss);
              field.nmiss   = (size_t)nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
380
381
	      field.grid    = vars1[varID].grid;
	      field.missval = vars1[varID].missval;

382
	      if ( operatorID == VERTINT && IS_NOT_EQUAL(layer_thickness, 1.0) ) farcmul(&field, layer_thickness);
383
	      if ( lmean && IS_NOT_EQUAL(layer_weight, 1.0) ) farcmul(&field, layer_weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
384

Uwe Schulzweida's avatar
Uwe Schulzweida committed
385
386
387
388
	      if ( field.nmiss > 0 || samp1[varID].ptr )
		{
		  if ( samp1[varID].ptr == NULL )
		    {
389
		      samp1[varID].ptr = (double*) Malloc(gridsize*sizeof(double));
390
		      for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
392
393
			samp1[varID].ptr[i] = vars1[varID].nsamp;
		    }

394
		  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
		    if ( !DBL_IS_EQUAL(field.ptr[i], vars1[varID].missval) )
396
		      samp1[varID].ptr[i] += layer_weight;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
398
		}

399
	      if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
		{
401
402
403
404
405
406
407
408
409
410
                  if ( IS_NOT_EQUAL(layer_weight, 1.0) )
                    {
                      farsumqw(&vars2[varID], field, layer_weight);
                      farsumw(&vars1[varID], field, layer_weight);
                    }
                  else
                    {
                      farsumq(&vars2[varID], field);
                      farsum(&vars1[varID], field);
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
413
414
415
416
              else if ( lrange )
                {
                  farmin(&vars2[varID], field);
                  farmax(&vars1[varID], field);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
418
419
420
421
422
423
424
425
	      else
		{
		  farfun(&vars1[varID], field, operfunc);
		}
	    }
	}

      for ( varID = 0; varID < nvars; varID++ )
	{
426
427
	  if ( vars1[varID].nsamp )
	    {
428
	      if ( lmean )
429
430
		{
		  if ( samp1[varID].ptr == NULL )
431
		    farcdiv(&vars1[varID], (double)vars1[varID].nsamp);
432
433
434
		  else
		    fardiv(&vars1[varID], samp1[varID]);
		}
435
	      else if ( lvarstd )
436
437
438
		{
		  if ( samp1[varID].ptr == NULL )
		    {
439
		      if ( lstd )
440
			farcstd(&vars1[varID], vars2[varID], vars1[varID].nsamp, divisor);
441
		      else
442
			farcvar(&vars1[varID], vars2[varID], vars1[varID].nsamp, divisor);
443
444
445
		    }
		  else
		    {
446
		      if ( lstd )
447
			farstd(&vars1[varID], vars2[varID], samp1[varID], divisor);
448
		      else
449
			farvar(&vars1[varID], vars2[varID], samp1[varID], divisor);
450
451
		    }
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452
453
454
455
              else if ( lrange )
                {
                  farsub(&vars1[varID], vars2[varID]);
                }
456
457

	      streamDefRecord(streamID2, varID, 0);
458
	      streamWriteRecord(streamID2, vars1[varID].ptr, (int)vars1[varID].nmiss);
459
460
	      vars1[varID].nsamp = 0;
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
461
462
463
464
465
466
467
	}

      tsID++;
    }

  for ( varID = 0; varID < nvars; varID++ )
    {
468
469
470
      Free(vars1[varID].ptr);
      if ( samp1[varID].ptr ) Free(samp1[varID].ptr);
      if ( lvarstd ) Free(vars2[varID].ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
471
472
    }

473
474
475
  Free(vars1);
  Free(samp1);
  if ( lvarstd ) Free(vars2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476

477
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
478

Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
  if ( needWeights )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
    for ( int index = 0; index < nzaxis; ++index )
481
      if ( vert[index].numlevel > 1 )  Free(vert[index].weights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482

483
484
  Free(vert);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
486
487
  streamClose(streamID2);
  streamClose(streamID1);

488
489
  vlistDestroy(vlistID2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
490
491
  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
492
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
}