Importbinary.cc 15.5 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-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
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.
*/

Oliver Heidmann's avatar
Oliver Heidmann committed
18
19
#include <cassert>

Ralf Mueller's avatar
Ralf Mueller committed
20
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
21

22
#include "cdo_options.h"
23
#include "process_int.h"
24
#include "cdo_cdi_wrapper.h"
Oliver Heidmann's avatar
Oliver Heidmann committed
25
#include "printinfo.h"
26
27
28
29
30

extern "C"
{
#include "lib/gradsdes/gradsdes.h"
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31

32
33
static void
get_dim_vals(dsets_t *pfi, double *vals, int dimlen, int dim)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
34
{
35
  gadouble (*conv)(gadouble *, gadouble);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
38
  gadouble *cvals;
  int i;

39
  assert(dimlen == pfi->dnum[dim]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40

41
  if (pfi->linear[dim] == 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
    {
43
44
45
46
47
      for (i = 0; i < dimlen; ++i)
        {
          vals[i] = pfi->grvals[dim][i + 1];
          /* printf("%d %g\n", i, vals[i]); */
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
    }
49
  else if (pfi->linear[dim] == 1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52
    {
      /*
      for ( i = 0; i < 3; ++i )
53
        printf("%d %g %g\n", i, pfi->grvals[dim][i] , pfi->abvals[dim][i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
56
      */
      conv = pfi->gr2ab[dim];
      cvals = pfi->grvals[dim];
57
58
59
60
61
      for (i = 0; i < dimlen; ++i)
        {
          vals[i] = conv(cvals, i + 1);
          /* printf("%d %g\n", i, vals[i]); */
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
63
64
    }
}

65
66
static void
rev_vals(double *vals, int n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
69
{
  double dum;

70
  for (int i = 0; i < n / 2; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
    {
72
      dum = vals[i];
73
74
      vals[i] = vals[n - 1 - i];
      vals[n - 1 - i] = dum;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
    }
}

78
static bool
79
y_is_gauss(Varray<double> &gridyvals, int ysize)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
{
81
  bool lgauss = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
  int i;

84
  if (ysize > 2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
    {
86
87
      Varray<double> yvals(ysize);
      Varray<double> yw(ysize);
88
      gaussianLatitudes(yvals.data(), yw.data(), (size_t) ysize);
89
      for (i = 0; i < ysize; i++) yvals[i] = std::asin(yvals[i]) / M_PI * 180.0;
90

Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
      for (i = 0; i < ysize; i++)
92
        if (std::fabs(yvals[i] - gridyvals[i]) > ((yvals[0] - yvals[1]) / 500)) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93

94
      if (i == ysize) lgauss = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95

96
      // check S->N
97
      if (lgauss == false)
98
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
          for (i = 0; i < ysize; i++)
100
            if (std::fabs(yvals[i] - gridyvals[ysize - i - 1]) > ((yvals[0] - yvals[1]) / 500)) break;
101

102
          if (i == ysize) lgauss = true;
103
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
    }

106
  return lgauss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
}

109
110
static int
define_grid(dsets_t *pfi)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
{
112
113
  const int nx = pfi->dnum[0];
  const int ny = pfi->dnum[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114

115
116
  Varray<double> xvals(nx);
  Varray<double> yvals(ny);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117

118
119
  get_dim_vals(pfi, xvals.data(), nx, 0);
  get_dim_vals(pfi, yvals.data(), ny, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120

121
  if (pfi->yrflg) rev_vals(yvals.data(), ny);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122

123
  const bool lgauss = (pfi->linear[1] == 0) ? y_is_gauss(yvals, (size_t) ny) : false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124

125
  const int gridtype = lgauss ? GRID_GAUSSIAN : GRID_LONLAT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126

127
  const int gridID = gridCreate(gridtype, nx * ny);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
129
130
  gridDefXsize(gridID, nx);
  gridDefYsize(gridID, ny);

131
132
  gridDefXvals(gridID, xvals.data());
  gridDefYvals(gridID, yvals.data());
133

134
  return gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
}

137
138
static int
define_level(dsets_t *pfi, int nlev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
{
  int zaxisID = -1;
141
  int nz = pfi->dnum[2];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142

143
  if (nz)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
    {
145
      Varray<double> zvals(nz);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146

147
      get_dim_vals(pfi, zvals.data(), nz, 2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148

149
150
      if (nz == 1 && IS_EQUAL(zvals[0], 0))
        zaxisID = zaxisCreate(ZAXIS_SURFACE, nz);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
      else
152
153
        {
          if (nlev > 0 && nlev < nz) nz = nlev;
154
          if (pfi->zrflg) rev_vals(zvals.data(), nz);
155
156
          zaxisID = zaxisCreate(ZAXIS_GENERIC, nz);
        }
157
      zaxisDefLevels(zaxisID, zvals.data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
158
159
160
161
162
163
164
165
166
167
    }
  else
    {
      double level = 0;
      nz = 1;

      zaxisID = zaxisCreate(ZAXIS_SURFACE, nz);
      zaxisDefLevels(zaxisID, &level);
    }

168
  return zaxisID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
}

171
172
void *
Importbinary(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
{
174
  size_t i;
175
  size_t nmiss = 0, n_nan;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
  int ivar;
177
  int varID = -1, levelID, tsID;
178
  int datatype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
  dsets_t pfi;
180
181
  int tcur, told, fnum;
  int tmin = 0, tmax = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182
  char *ch = nullptr;
183
  int nlevels;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
  int e, flag;
185
  size_t rc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
  struct dt dtim, dtimi;
187
  double missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
  double fmin, fmax;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
  double sfclevel = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190

191
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192

193
194
  operatorCheckArgc(0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
  dsets_init(&pfi);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
196

197
  int status = read_gradsdes((char *) cdoGetStreamName(0), &pfi);
198
  if (Options::cdoVerbose) fprintf(stderr, "status %d\n", status);
199
200
  // if ( status ) cdoAbort("Open failed on %s!", pfi.name);
  if (status) cdoAbort("Open failed!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201

202
203
  int nrecs = pfi.trecs;
  int nvars = pfi.vnum;
204
  struct gavar *pvar = pfi.pvar1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205

206
  if (nvars == 0) cdoAbort("No variables found!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207

208
209
  const auto gridID = define_grid(&pfi);
  const auto zaxisID = define_level(&pfi, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210

211
  const auto zaxisIDsfc = zaxisCreate(ZAXIS_SURFACE, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
  zaxisDefLevels(zaxisIDsfc, &sfclevel);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213

214
  const auto vlistID = vlistCreate();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215

216
217
  Varray<int> var_zaxisID(nvars);
  Varray<int> var_dfrm(nrecs);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
  std::vector<RecordInfo> recList(nrecs);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219

220
  int recID = 0;
221
  for (ivar = 0; ivar < nvars; ++ivar)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
223
    {
      /*
224
      if ( Options::cdoVerbose )
225
226
227
        fprintf(stderr, "1:%s 2:%s %d %d %d %d 3:%s %d \n",
                pvar->abbrv, pvar->longnm, pvar->offset, pvar->recoff,
      pvar->levels, pvar->nvardims, pvar->varnm, pvar->var_t);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
      nlevels = pvar->levels;
230
231
232
233
234
235

      if (nlevels == 0)
        {
          nlevels = 1;
          varID = vlistDefVar(vlistID, gridID, zaxisIDsfc, TIME_VARYING);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
      else
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
        {
          if (nlevels > zaxisInqSize(zaxisID))
            cdoAbort("Variable %s has too many number of levels!", pvar->abbrv);
          else if (nlevels < zaxisInqSize(zaxisID))
            {
              int vid, zid = -1, nlev;
              for (vid = 0; vid < ivar; ++vid)
                {
                  zid = var_zaxisID[vid];
                  nlev = zaxisInqSize(zid);
                  if (nlev == nlevels) break;
                }

              if (vid == ivar) zid = define_level(&pfi, nlevels);
              varID = vlistDefVar(vlistID, gridID, zid, TIME_VARYING);
            }
          else
            varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
256

257
258
      var_zaxisID[varID] = vlistInqVarZaxis(vlistID, varID);

259
      cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, pvar->abbrv);
260
      {
261
        char *longname = pvar->varnm;
262
        const auto len = strlen(longname);
263
264
265
266
267
268
        if (longname[0] == '\'' && longname[len - 1] == '\'')
          {
            longname[len - 1] = 0;
            longname++;
          }
        if (longname[0] == '\t') longname++;
269
        cdiDefKeyString(vlistID, varID, CDI_KEY_LONGNAME, longname);
270
      }
271

272
      missval = pfi.undef;
273
      datatype = CDI_DATATYPE_FLT32;
274

275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
      if (pvar->dfrm == 1)
        {
          datatype = CDI_DATATYPE_UINT8;
          if (missval < 0 || missval > 255) missval = 255;
        }
      else if (pvar->dfrm == 2)
        {
          datatype = CDI_DATATYPE_UINT16;
          if (missval < 0 || missval > 65535) missval = 65535;
        }
      else if (pvar->dfrm == -2)
        {
          datatype = CDI_DATATYPE_INT16;
          if (missval < -32768 || missval > 32767) missval = -32768;
        }
      else if (pvar->dfrm == 4)
        {
          datatype = CDI_DATATYPE_INT32;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293
          if (missval < -2147483648 || missval > 2147483647) missval = -2147483646;
294
295
296
297
        }
      else if (pfi.flt64)
        datatype = CDI_DATATYPE_FLT64;

298
299
      vlistDefVarDatatype(vlistID, varID, datatype);
      vlistDefVarMissval(vlistID, varID, missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
300

301
302
      for (levelID = 0; levelID < nlevels; ++levelID)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
303
          if (recID >= nrecs) cdoAbort("Internal problem with number of records!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
305
          recList[recID].varID = varID;
          recList[recID].levelID = levelID;
306
307
308
          var_dfrm[recID] = pvar->dfrm;
          recID++;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312

      pvar++;
    }

313
314
  const auto calendar = CALENDAR_STANDARD;
  const auto taxisID = cdoTaxisCreate(TAXIS_RELATIVE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
  taxisDefCalendar(taxisID, calendar);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
317
  vlistDefTaxis(vlistID, taxisID);

318
  const auto streamID = cdoOpenWrite(1);
319
  cdoDefVlist(streamID, vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320

321
322
  const size_t gridsize = pfi.dnum[0] * pfi.dnum[1];
  int recoffset = pfi.xyhdr * (pfi.flt64 ? 8 : 4);
323
  if (pfi.seqflg) recoffset += 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324

325
  // recsize = pfi.gsiz*4;
326
327
  size_t recsize = pfi.gsiz * 8;
  Varray<char> rec(recsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328

329
  Varray<double> array(gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330

Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
332
333
334
  if (pfi.tmplat)
    for ( i = 0; i <  pfi.dnum[3]; ++i )
      printf("%d %d\n", i, pfi.fnums[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336

Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
  pfi.infile = nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
  tcur = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
  e = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
  while (1)
341
    { /* loop over all times for this ensemble */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
      if (pfi.tmplat)
343
344
        {
          /* make sure no file is open */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
345
          if (pfi.infile != nullptr)
346
347
            {
              fclose(pfi.infile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
              pfi.infile = nullptr;
349
350
351
352
353
354
            }
          /* advance to first valid time step for this ensemble */
          if (tcur == 0)
            {
              told = 0;
              tcur = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
355
              while (pfi.fnums[tcur - 1] == -1) tcur++;
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
            }
          else
            { /* tcur!=0 */
              told = pfi.fnums[tcur - 1];
              /* increment time step until fnums changes */
              while (told == pfi.fnums[tcur - 1] && tcur <= pfi.dnum[3])
                {
                  tcur++;
                  if (tcur > pfi.dnum[3]) break;
                }
            }

          /* make sure we haven't advanced past end of time axis */
          if (tcur > pfi.dnum[3]) break;

          /* check if we're past all valid time steps for this ensemble */
          if ((told != -1) && (pfi.fnums[tcur - 1] == -1)) break;

          /* Find the range of t indexes that have the same fnums value.
             These are the times that are contained in this particular file */
          tmin = tcur;
          tmax = tcur - 1;
          fnum = pfi.fnums[tcur - 1];
          if (fnum != -1)
            {
              while (fnum == pfi.fnums[tmax])
                {
                  tmax++;
                  if (tmax == pfi.dnum[3]) break;
                }
              gr2t(pfi.grvals[3], (gadouble) tcur, &dtim);
              gr2t(pfi.grvals[3], (gadouble) 1, &dtimi);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
388
389
              ch = gafndt(pfi.name, &dtim, &dtimi, pfi.abvals[3], pfi.pchsub1, nullptr, tcur, e, &flag);
              if (ch == nullptr) cdoAbort("Couldn't determine data file name for e=%d t=%d!", e, tcur);
390
391
392
393
394
395
396
397
398
399
            }
        }
      else
        {
          /* Data set is not templated */
          ch = pfi.name;
          tmin = 1;
          tmax = pfi.dnum[3];
        }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
      /* Open this file and position to start of first record */
401
      if (Options::cdoVerbose) cdoPrint("Opening file: %s", ch);
402
      pfi.infile = fopen(ch, "rb");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403
      if (pfi.infile == nullptr)
404
405
406
        {
          if (pfi.tmplat)
            {
407
              cdoWarning("Could not open file: %s", ch);
408
409
410
411
412
413
414
              break;
            }
          else
            {
              cdoAbort("Could not open file: %s", ch);
            }
        }
415
416
417

      /* file header */
      if (pfi.fhdr > 0) fseeko(pfi.infile, pfi.fhdr, SEEK_SET);
418

Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
      /* Get file size */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
      /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
422
      fseeko(pfi.infile,0L,2);
      flen = ftello(pfi.infile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
423
424

      printf("flen %d tsiz %d\n", flen, pfi.tsiz);
425

Uwe Schulzweida's avatar
Uwe Schulzweida committed
426
      fseeko (pfi.infile,0,0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
      */
428
429
430
      for (tsID = tmin - 1; tsID < tmax; ++tsID)
        {
          gr2t(pfi.grvals[3], (gadouble)(tsID + 1), &dtim);
431
432
          const auto vdate = cdiEncodeDate(dtim.yr, dtim.mo, dtim.dy);
          const auto vtime = cdiEncodeTime(dtim.hr, dtim.mn, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
434
          if (Options::cdoVerbose)
            cdoPrint(" Reading timestep: %3d %s %s", tsID + 1, dateToString(vdate).c_str(), timeToString(vtime).c_str());
435
436
437

          taxisDefVdate(taxisID, vdate);
          taxisDefVtime(taxisID, vtime);
438
          cdoDefTimestep(streamID, tsID);
439

440
          for (recID = 0; recID < nrecs; ++recID)
441
442
443
444
445
446
447
448
449
450
451
452
            {
              /* record size depends on data type */
              if (var_dfrm[recID] == 1)
                {
                  recsize = pfi.gsiz;
                }
              else if ((var_dfrm[recID] == 2) || (var_dfrm[recID] == -2))
                {
                  recsize = pfi.gsiz * 2;
                }
              else
                {
453
                  recsize = pfi.flt64 ? pfi.gsiz * 8 : pfi.gsiz * 4;
454
455
                }

456
              rc = fread(rec.data(), 1, recsize, pfi.infile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
              if (rc < recsize) cdoAbort("I/O error reading record=%d of timestep=%d!", recID + 1, tsID + 1);
458
459
460
461

              /* convert */
              if (var_dfrm[recID] == 1)
                {
462
                  unsigned char *carray = (unsigned char *) &rec[recoffset];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
463
                  for (i = 0; i < gridsize; ++i) array[i] = (double) carray[i];
464
465
466
                }
              else if (var_dfrm[recID] == 2)
                {
467
                  unsigned short *sarray = (unsigned short *) &rec[recoffset];
468
                  if (pfi.bswap) gabswp2(sarray, gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
                  for (i = 0; i < gridsize; ++i) array[i] = (double) sarray[i];
470
471
472
                }
              else if (var_dfrm[recID] == -2)
                {
473
                  short *sarray = (short *) &rec[recoffset];
474
                  if (pfi.bswap) gabswp2(sarray, gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
                  for (i = 0; i < gridsize; ++i) array[i] = (double) sarray[i];
476
477
478
                }
              else if (var_dfrm[recID] == 4)
                {
479
                  int *iarray = (int *) &rec[recoffset];
480
                  if (pfi.bswap) gabswp(iarray, gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
                  for (i = 0; i < gridsize; ++i) array[i] = (double) iarray[i];
482
483
484
485
486
                }
              else
                {
                  if (pfi.flt64)
                    {
487
                      double *darray = (double *) &rec[recoffset];
488
                      if (pfi.bswap) gabswp(darray, gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
                      for (i = 0; i < gridsize; ++i) array[i] = darray[i];
490
491
492
                    }
                  else
                    {
493
                      float *farray = (float *) &rec[recoffset];
494
                      if (pfi.bswap) gabswp(farray, gridsize);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
                      for (i = 0; i < gridsize; ++i) array[i] = (double) farray[i];
496
497
498
499
500
501
502
503
504
505
506
507
508
509
                    }
                }

              fmin = 1.e99;
              fmax = -1.e99;
              nmiss = 0;
              n_nan = 0;
              for (i = 0; i < gridsize; ++i)
                {
                  if (array[i] > pfi.ulow && array[i] < pfi.uhi)
                    {
                      array[i] = pfi.undef;
                      nmiss++;
                    }
510
                  else if (std::isnan(array[i]))
511
512
513
514
515
516
517
                    {
                      array[i] = pfi.undef;
                      nmiss++;
                      n_nan++;
                    }
                  else
                    {
518
519
                      fmin = std::min(fmin, array[i]);
                      fmax = std::max(fmax, array[i]);
520
521
522
                    }
                }
              /*
523
              if ( Options::cdoVerbose )
524
525
526
                printf("%3d %4d %3d %6zu %6zu %12.5g %12.5g\n", tsID, recID,
              recoffset, nmiss, n_nan, fmin, fmax);
              */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
528
              varID = recList[recID].varID;
              levelID = recList[recID].levelID;
529
              cdoDefRecord(streamID, varID, levelID);
530
              cdoWriteRecord(streamID, array.data(), nmiss);
531
532
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533

Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
535
      /* break out if not templating */
      if (!pfi.tmplat) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536

537
    } /* end of while (1) loop */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
538

539
  processDefVarNum(vlistNvars(vlistID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
540

541
  cdoStreamClose(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
543
544
545
546
547

  vlistDestroy(vlistID);
  gridDestroy(gridID);
  zaxisDestroy(zaxisID);
  taxisDestroy(taxisID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
  if (pfi.infile) fclose(pfi.infile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
551

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
553
}