Gradsdes.cc 38.2 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-2018 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:

21
      Gradsdes   gradsdes       GrADS data descriptor file
22
23
*/

Ralf Mueller's avatar
Ralf Mueller committed
24
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
25

Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
#include "cdo_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
#include "grid.h"
28
#include "pstream_int.h"
29
#include "util_fileextensions.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30

31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
/*
  Values output into the grib map file:

  Header:

  hipnt info:  0 - version number (1)
               1 - number of times in file
               2 - number of records per time
               3 - Grid type
                 255 - user defined grid.  descriptor
                       describes grid exactly; one record
                       per grid.
                  29 - Predefined grid set 29 and 30.
                       Two records per grid.

  hfpnt info:  None

  Info:

  intpnt info (for each mapped grib record) :
                 0 - position of start of data in file
                 1 - position of start of bit map in file
                 2 - number of bits per data element

  fltpnt info :
                 0 - decimal scale factor for this record
                 1 - binary scale factor
                 2 - reference value

*/
61
62
63
64
65
66
67
68
69
70
71
struct gaindx
{
  int type;      /* Indexing file type             */
  int hinum;     /* Number of ints in header       */
  int hfnum;     /* Number of floats in header     */
  int intnum;    /* Number of index ints (long)    */
  int fltnum;    /* Number of index floats         */
  int *hipnt;    /* Pointer to header int values   */
  float *hfpnt;  /* Pointer to header float values */
  int *intpnt;   /* Pointer to int index values    */
  float *fltpnt; /* Pointer to float index values  */
72
};
73
74
75
76
struct gaindxb
{
  int bignum;    /* Number of off_t values */
  off_t *bigpnt; /* Pointer to off_t values */
77
};
78

Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80
/* Byte swap requested number of 4 byte elements */

81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
static void
gabswp(void *r, int cnt)
{
  int i;
  char *ch1, *ch2, *ch3, *ch4, cc1, cc2;

  ch1 = (char *) r;
  ch2 = ch1 + 1;
  ch3 = ch2 + 1;
  ch4 = ch3 + 1;
  for (i = 0; i < cnt; i++)
    {
      cc1 = *ch1;
      cc2 = *ch2;
      *ch1 = *ch4;
      *ch2 = *ch3;
      *ch3 = cc2;
      *ch4 = cc1;
      ch1 += 4;
      ch2 += 4;
      ch3 += 4;
      ch4 += 4;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
106
}

/*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
 * convert an IBM float to single precision number v1.0
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
 *
 *                      Wesley Ebisuzaki
 */
111
112
113
static float
ibm2flt(unsigned char *ibm)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114

115
116
117
118
119
120
121
122
123
124
125
126
127
  int positive, power;
  unsigned int abspower;
  long int mant;
  double value, exp;

  positive = (ibm[0] & 0x80) == 0;
  mant = (ibm[1] << 16) + (ibm[2] << 8) + ibm[3];
  power = (int) (ibm[0] & 0x7f) - 64;
  abspower = power > 0 ? power : -power;

  /* calc exp */
  exp = 16.0;
  value = 1.0;
128
129
130
131
132
133
134
135
  while (abspower)
    {
      if (abspower & 1)
        {
          value *= exp;
        }
      exp = exp * exp;
      abspower >>= 1;
136
137
138
139
140
    }

  if (power < 0) value = 1.0 / value;
  value = value * mant / 16777216.0;
  if (positive == 0) value = -value;
141
  return (float) value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
144
145
146
147
148
149
150
}

/*
 * convert a float to an IBM single precision number v1.0
 *
 *                      Wesley Ebisuzaki
 *
 * doesn't handle subnormal numbers
 */
151
152
153
static int
flt2ibm(float x, unsigned char *ibm)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154

155
156
157
  int sign, exp, i;
  double mant;

158
159
160
161
162
  if (!(fabs((double) x) > 0))
    {
      ibm[0] = ibm[1] = ibm[2] = ibm[3] = 0;
      return 0;
    }
163
164

  /* sign bit */
165
166
167
168
169
170
171
  if (x < 0.0)
    {
      sign = 128;
      x = -x;
    }
  else
    sign = 0;
172
173
174
175
176
177

  mant = frexp((double) x, &exp);

  /* round up by adding 2**-24 */
  /* mant = mant + 1.0/16777216.0; */

178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  if (mant >= 1.0)
    {
      mant = 0.5;
      exp++;
    }
  while (exp & 3)
    {
      mant *= 0.5;
      exp++;
    }

  exp = exp / 4 + 64;

  if (exp < 0)
    {
      fprintf(stderr, "underflow in flt2ibm\n");
      ibm[0] = ibm[1] = ibm[2] = ibm[3] = 0;
      return 0;
    }
  if (exp > 127)
    {
      fprintf(stderr, "overflow in flt2ibm\n");
      ibm[0] = sign | 127;
      ibm[1] = ibm[2] = ibm[3] = 255;
      return -1;
    }
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218

  /* normal number */

  ibm[0] = sign | exp;

  mant = mant * 256.0;
  i = (int) floor(mant);
  mant = mant - i;
  ibm[1] = i;

  mant = mant * 256.0;
  i = (int) floor(mant);
  mant = mant - i;
  ibm[2] = i;

219
  ibm[3] = (int) floor(mant * 256.0);
220
221

  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
223
}

224
225
#define GET_UINT4(a, b, c, d) ((int) ((a << 24) + (b << 16) + (c << 8) + (d)))
#define Put1Byte(buf, cnt, ival) (buf[cnt++] = (ival))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
228
#define Put2Byte(buf, cnt, ival) ((buf[cnt++] = (ival) >> 8), (buf[cnt++] = (ival)))
#define Put4Byte(buf, cnt, ival) \
  ((buf[cnt++] = (ival) >> 24), (buf[cnt++] = (ival) >> 16), (buf[cnt++] = (ival) >> 8), (buf[cnt++] = (ival)))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229

Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
#define PutInt(buf, cnt, ival) (ival < 0 ? Put4Byte(buf, cnt, 0x7fffffff - ival + 1) : Put4Byte(buf, cnt, ival))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231

232
233
static void
dumpmap()
234
235
236
237
238
239
{
  unsigned char urec[4];
  unsigned char vermap;
  unsigned char mrec[512];
  int swpflg = 0;
  int i;
240
  int nrecords = 0;
241
  struct gaindx indx;
242
  struct gaindxb indxb;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
  size_t nbytes;
244
245
246
247
248
249
  FILE *mapfp;

  indx.hipnt = NULL;
  indx.hfpnt = NULL;
  indx.intpnt = NULL;
  indx.fltpnt = NULL;
250
251
  indxb.bigpnt = NULL;
  indxb.bignum = 0;
252

253
254
  mapfp = fopen(cdoGetStreamName(0), "r");
  if (mapfp == NULL) cdoAbort("Open failed on %s", cdoGetStreamName(0));
255
256
257
258

  /* check the version number */

  fseek(mapfp, 1, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
  nbytes = fread(&vermap, sizeof(unsigned char), 1, mapfp);
260

261
  if (vermap == 0) vermap = 1;
262
263
264

  printf("gribmap version = %d\n", vermap);

265
  if (vermap == 2)
266
267
268
    {
      fseek(mapfp, 2, 0);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
270
      indx.hinum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
271

Uwe Schulzweida's avatar
Uwe Schulzweida committed
272
      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
273
      indx.hfnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
274

Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
276
      indx.intnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
277

Uwe Schulzweida's avatar
Uwe Schulzweida committed
278
      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
279
      indx.fltnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
280

Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
      nbytes = fread(mrec, sizeof(unsigned char), 7, mapfp);
282

283
      if (indx.hinum > 0)
284
        {
285
286
          indx.hipnt = (int *) Malloc(sizeof(int) * indx.hinum);
          for (i = 0; i < indx.hinum; i++)
287
288
            {
              nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
289
              indx.hipnt[i] = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
290
291
            }
        }
292
      if (indx.hfnum > 0)
293
        {
294
295
          indx.hfpnt = (float *) Malloc(sizeof(float) * indx.hfnum);
          nbytes = fread(indx.hfpnt, sizeof(float), indx.hfnum, mapfp);
296
        }
297
      if (indx.intnum > 0)
298
        {
299
300
          indx.intpnt = (int *) Malloc(sizeof(int) * indx.intnum);
          for (i = 0; i < indx.intnum; i++)
301
302
            {
              nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
303
              indx.intpnt[i] = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
              if (indx.intpnt[i] < 0) indx.intpnt[i] = 0x7fffffff - indx.intpnt[i] + 1;
305
306
            }
        }
307
      if (indx.fltnum > 0)
308
        {
309
310
          indx.fltpnt = (float *) Malloc(sizeof(float) * indx.fltnum);
          for (i = 0; i < indx.fltnum; i++)
311
312
313
314
315
            {
              nbytes = fread(urec, sizeof(unsigned char), 4, mapfp);
              indx.fltpnt[i] = ibm2flt(urec);
            }
        }
316
317
318
319
    }
  else
    {
      fseek(mapfp, 0, 0);
320
      nbytes = fread(&indx, sizeof(struct gaindx), 1, mapfp);
321

322
323
324
      if (indx.type >> 24 > 0) swpflg = 1;
      if (swpflg) printf("swap endian!\n");
      if (swpflg) gabswp((float *) &indx.type, 5);
325

326
      if (indx.hinum > 0)
327
        {
328
329
330
          indx.hipnt = (int *) Malloc(sizeof(int) * indx.hinum);
          nbytes = fread(indx.hipnt, sizeof(int), indx.hinum, mapfp);
          if (swpflg) gabswp((float *) (indx.hipnt), indx.hinum);
331
        }
332
      if (indx.hfnum > 0)
333
        {
334
335
336
          indx.hfpnt = (float *) Malloc(sizeof(float) * indx.hfnum);
          nbytes = fread(indx.hfpnt, sizeof(float), indx.hfnum, mapfp);
          if (swpflg) gabswp(indx.hfpnt, indx.hfnum);
337
        }
338

339
      if (indx.intnum > 0)
340
        {
341
342
343
          indx.intpnt = (int *) Malloc(sizeof(int) * indx.intnum);
          nbytes = fread(indx.intpnt, sizeof(int), indx.intnum, mapfp);
          if (swpflg) gabswp((float *) (indx.intpnt), indx.intnum);
344
        }
345
      if (indx.fltnum > 0)
346
        {
347
348
349
          indx.fltpnt = (float *) Malloc(sizeof(float) * indx.fltnum);
          nbytes = fread(indx.fltpnt, sizeof(float), indx.fltnum, mapfp);
          if (swpflg) gabswp(indx.fltpnt, indx.fltnum);
350
        }
351

352
      if (indx.hipnt[0] == 4)
353
354
        {
          indxb.bignum = indx.hipnt[4];
355
          if (indxb.bignum > 0)
356
            {
357
358
359
              indxb.bigpnt = (off_t *) Malloc(sizeof(off_t) * indxb.bignum);
              nbytes = fread(indxb.bigpnt, sizeof(off_t), indxb.bignum, mapfp);
              if (swpflg) gabswp(indxb.bigpnt, indxb.bignum);
360
361
            }
        }
362
363
    }

364
365
  UNUSED(nbytes);

366
367
368
  fclose(mapfp);

  printf("hinum: %d\n", indx.hinum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
  for (i = 0; i < indx.hinum; i++) printf("%3d %5d\n", i + 1, indx.hipnt[i]);
370

371
372
  printf("\n");
  printf("hfnum: %d\n", indx.hfnum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
  for (i = 0; i < indx.hfnum; i++) printf("%3d %g\n", i + 1, indx.hfpnt[i]);
374

375
  printf("\n");
376

377
  nrecords = indx.hipnt[1] * indx.hipnt[2];
378

379
  if (indx.intnum == indx.fltnum)
380
381
    {
      printf("num: %d\n", indx.intnum);
382
      for (i = 0; i < indx.intnum / 3; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
384
        printf("%3d %8d %6d %4d %8g %10g %8g\n", i + 1, indx.intpnt[i * 3], indx.intpnt[i * 3 + 1], indx.intpnt[i * 3 + 2],
               indx.fltpnt[i * 3], indx.fltpnt[i * 3 + 1], indx.fltpnt[i * 3 + 2]);
385
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
  else if (indx.intnum == nrecords && indx.fltnum == nrecords * 3 && indxb.bignum == nrecords * 2)
387
388
    {
      printf("nrecords: %d\n", nrecords);
389
      for (i = 0; i < nrecords; i++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
391
        printf("%3d %8zd %6zd %4d %8g %10g %8g\n", i + 1, (size_t) indxb.bigpnt[i * 2], (size_t) indxb.bigpnt[i * 2 + 1],
               indx.intpnt[i], indx.fltpnt[i * 3], indx.fltpnt[i * 3 + 1], indx.fltpnt[i * 3 + 2]);
392
    }
393
394
395
  else
    {
      printf("intnum: %d\n", indx.intnum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
      for (i = 0; i < indx.intnum; i++) printf("%3d %d\n", i + 1, indx.intpnt[i]);
397
398
399

      printf("\n");
      printf("fltnum: %d\n", indx.fltnum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
      for (i = 0; i < indx.fltnum; i++) printf("%3d %g\n", i + 1, indx.fltpnt[i]);
401
402
403

      printf("\n");
      printf("bignum: %d\n", indxb.bignum);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404
      for (i = 0; i < indxb.bignum; i++) printf("%3d %zd\n", i + 1, (size_t) indxb.bigpnt[i]);
405
406
407
    }
}

408
409
static void
ctl_xydef(FILE *gdp, int gridID, bool *yrev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
411
412
413
{
  int i, j;
  double xfirst, yfirst, xinc, yinc;

414
  *yrev = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415

416
417
  int xsize = gridInqXsize(gridID);
  int ysize = gridInqYsize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418

419
420
  int gridtype = gridInqType(gridID);
  int projtype = gridInqProjType(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
422
423

  /* XDEF */

424
  if (gridtype == GRID_PROJECTION && projtype == CDI_PROJ_LCC)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
426
427
428
    {
      double xmin = 1.e10, xmax = -1.e10, ymin = 1.e10, ymax = -1.e10;
      double xrange, yrange;
      int nx, ny, ni;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
429
      double inc[] = { 1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001 };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430

431
432
433
      double xinc = gridInqXinc(gridID);
      double yinc = gridInqYinc(gridID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
434
      double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;
435

Uwe Schulzweida's avatar
Uwe Schulzweida committed
436
      gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437
      fprintf(gdp, "PDEF %d %d LCCR %g %g 1 1 %g %g %g %g %g\n", xsize, ysize, xval_0, yval_0, lat_1, lat_2, lon_0, xinc, yinc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
439

      gridID = gridToCurvilinear(gridID, 0);
440
441
442
443
      std::vector<double> xvals(xsize * ysize);
      std::vector<double> yvals(xsize * ysize);
      gridInqXvals(gridID, xvals.data());
      gridInqYvals(gridID, yvals.data());
444
      for (i = 0; i < xsize * ysize; ++i)
445
        {
446
447
448
449
450
          if (xvals[i] > 180) xvals[i] -= 360;
          if (xvals[i] < xmin) xmin = xvals[i];
          if (xvals[i] > xmax) xmax = xvals[i];
          if (yvals[i] < ymin) ymin = yvals[i];
          if (yvals[i] > ymax) ymax = yvals[i];
451
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452

453
454
455
456
      xfirst = ((int) (xmin - 0.0));
      yfirst = ((int) (ymin - 0.0));
      xrange = ((int) (xmax + 1.5)) - xfirst;
      yrange = ((int) (ymax + 1.5)) - yfirst;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457

458
459
      ni = sizeof(inc) / sizeof(inc[0]);
      for (i = 0; i < ni; i++)
460
461
462
463
        {
          xinc = yinc = inc[i];
          nx = 1 + (int) (xrange / xinc);
          ny = 1 + (int) (yrange / yinc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464

465
          if (nx > 1.5 * xsize && ny > 1.5 * ysize) break;
466
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
468
469
470
471
472
473
474
475
476

      fprintf(gdp, "XDEF %d LINEAR %f %f\n", nx, xfirst, xinc);
      fprintf(gdp, "YDEF %d LINEAR %f %f\n", ny, yfirst, yinc);

      fprintf(gdp, "* XDEF 3600 LINEAR -179.95 0.1\n");
      fprintf(gdp, "* YDEF 1800 LINEAR  -89.95 0.1\n");
    }
  else
    {
      xfirst = gridInqXval(gridID, 0);
477
478
      xinc = gridInqXinc(gridID);
      if (IS_EQUAL(xinc, 0) && gridInqXvals(gridID, NULL))
479
        {
480
481
          std::vector<double> xvals(xsize);
          gridInqXvals(gridID, xvals.data());
482
          fprintf(gdp, "XDEF %d LEVELS ", xsize);
483
          j = 0;
484
          for (i = 0; i < xsize; i++)
485
486
487
            {
              fprintf(gdp, "%7.3f ", xvals[i]);
              j++;
488
              if (j == 6)
489
490
491
                {
                  fprintf(gdp, "\n");
                  j = 0;
492
                  if (i != xsize - 1) fprintf(gdp, "               ");
493
494
                }
            }
495
          if (j) fprintf(gdp, "\n");
496
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
      else
498
        {
499
          if (IS_EQUAL(xinc, 0)) xinc = 360.0 / xsize;
500
501
          fprintf(gdp, "XDEF %d LINEAR %f %f\n", xsize, xfirst, xinc);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
503
504
505
    }

  /* YDEF */

506
  if (!(gridtype == GRID_PROJECTION && projtype == CDI_PROJ_LCC))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
508
    {
      yfirst = gridInqYval(gridID, 0);
509
510
      yinc = gridInqYinc(gridID);
      if (gridtype == GRID_GAUSSIAN) yinc = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511

512
      if (IS_EQUAL(yinc, 0) && gridInqYvals(gridID, NULL))
513
        {
514
515
          std::vector<double> yvals(ysize);
          gridInqYvals(gridID, yvals.data());
516
          fprintf(gdp, "YDEF %d LEVELS ", ysize);
517
          j = 0;
518
          if (yvals[0] > yvals[ysize - 1])
519
            {
520
              *yrev = true;
521
              for (i = ysize - 1; i >= 0; i--)
522
523
524
                {
                  fprintf(gdp, "%7.3f ", yvals[i]);
                  j++;
525
                  if (j == 6)
526
527
528
                    {
                      fprintf(gdp, "\n");
                      j = 0;
529
                      if (i != 0) fprintf(gdp, "               ");
530
531
532
533
534
                    }
                }
            }
          else
            {
535
              for (i = 0; i < ysize; i++)
536
537
538
                {
                  fprintf(gdp, "%7.3f ", yvals[i]);
                  j++;
539
                  if (j == 6)
540
541
542
                    {
                      fprintf(gdp, "\n");
                      j = 0;
543
                      if (i != ysize - 1) fprintf(gdp, "               ");
544
545
546
547
                    }
                }
            }

548
          if (j) fprintf(gdp, "\n");
549
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
550
      else
551
        {
552
553
          if (IS_EQUAL(yinc, 0)) yinc = 180.0 / ysize;
          if (yinc < 0)
554
            {
555
              *yrev = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
556
              fprintf(gdp, "YDEF %d LINEAR %f %f\n", ysize, yfirst + yinc * (ysize - 1), -yinc);
557
558
559
560
            }
          else
            fprintf(gdp, "YDEF %d LINEAR %f %f\n", ysize, yfirst, yinc);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561
562
563
    }
}

564
565
static void
ctl_zdef(FILE *gdp, int vlistID, bool *zrev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
566
567
{
  int i, j, index;
568
569
570
  int zaxisIDmax = -1;
  int zaxisID, nlev;
  double levinc = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571

572
  *zrev = false;
573
  int nzaxis = vlistNzaxis(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574

575
  int nlevmax = 0;
576
  for (index = 0; index < nzaxis; index++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
    {
      zaxisID = vlistZaxis(vlistID, index);
579
580
      nlev = zaxisInqSize(zaxisID);
      if (nlev > nlevmax)
581
582
583
584
        {
          nlevmax = nlev;
          zaxisIDmax = zaxisID;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
585
586
    }

587
588
  std::vector<double> levels(nlevmax);
  cdoZaxisInqLevels(zaxisIDmax, levels.data());
589
590
  bool lplev = (zaxisInqType(zaxisIDmax) == ZAXIS_PRESSURE);
  double level0 = levels[0];
591
  if (nlevmax > 1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
      if (levels[0] < levels[1] && zaxisInqType(zaxisIDmax) != ZAXIS_HYBRID) *zrev = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
594
595
596

      levinc = levels[1] - levels[0];

597
      if (IS_EQUAL(levinc, 1)) *zrev = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598

599
      for (i = 1; i < nlevmax; i++)
600
        {
601
          if (IS_NOT_EQUAL(levinc, (levels[i] - levels[i - 1])))
602
603
604
605
606
            {
              levinc = 0;
              break;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
    }

609
610
  if (IS_NOT_EQUAL(levinc, 0))
    fprintf(gdp, "ZDEF %d LINEAR %g %g\n", nlevmax, level0, levinc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
612
613
  else
    {
      fprintf(gdp, "ZDEF %d LEVELS ", nlevmax);
614
      j = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
616
      /* zrev not needed !!!
      if ( *zrev )
617
618
619
620
621
622
623
624
625
626
627
628
629
630
        {
          for ( i = nlevmax-1; i >=0 ; i-- )
            {
              if ( lplev ) fprintf(gdp, "%g ", levels[i]/100);
              else         fprintf(gdp, "%d ", (int) levels[i]);
              j++;
              if ( j == 10 )
                {
                  fprintf(gdp, "\n");
                  j = 0;
                  if ( i != 0 ) fprintf(gdp, "               ");
                }
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631
632
      else
      */
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
      {
        for (i = 0; i < nlevmax; i++)
          {
            if (lplev)
              fprintf(gdp, "%g ", levels[i] / 100);
            else
              fprintf(gdp, "%g ", levels[i]);
            j++;
            if (j == 10)
              {
                fprintf(gdp, "\n");
                j = 0;
                if (i != (nlevmax - 1)) fprintf(gdp, "               ");
              }
          }
      }
      if (j) fprintf(gdp, "\n");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
650
651
652
    }
}

653
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
ctl_options(FILE *gdp, bool yrev, bool zrev, bool sequential, bool bigendian, bool littleendian, bool flt64, bool cal365day)
655
{
656
  /* if ( filetype == CDI_FILETYPE_GRB ) zrev = false; */
657

658
  if (yrev || zrev || sequential || bigendian || littleendian || flt64)
659
660
    {
      fprintf(gdp, "OPTIONS");
661
662
663
664
665
666
667
      if (yrev) fprintf(gdp, " yrev");
      if (zrev) fprintf(gdp, " zrev");
      if (sequential) fprintf(gdp, " sequential");
      if (bigendian) fprintf(gdp, " big_endian");
      if (littleendian) fprintf(gdp, " little_endian");
      if (flt64) fprintf(gdp, " flt64");
      if (cal365day) fprintf(gdp, " 365_day_calendar");
668
669
670
671
      fprintf(gdp, "\n");
    }
}

672
673
static void
ctl_undef(FILE *gdp, int vlistID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
674
{
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
675
  double missval = vlistInqVarMissval(vlistID, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
677
678
  fprintf(gdp, "UNDEF  %g\n", missval);
}

679
680
static void
ctl_vars(FILE *gdp, int filetype, int vlistID, int nvarsout, int *vars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
681
682
683
684
685
686
687
{
  int ltype, code;
  int zaxisID, nlev;
  int i, j;
  int len;
  char varname[CDI_MAX_NAME], varlongname[CDI_MAX_NAME], varunits[CDI_MAX_NAME];

Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
688
  int nvars = vlistNvars(vlistID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
689
690
691

  fprintf(gdp, "VARS  %d\n", nvarsout);

692
  for (int varID = 0; varID < nvars; varID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693
    {
694
      if (vars[varID] == TRUE)
695
696
        {
          zaxisID = vlistInqVarZaxis(vlistID, varID);
697
698
          ltype = zaxisInqLtype(zaxisID);
          nlev = zaxisInqSize(zaxisID);
699
700
701
          vlistInqVarName(vlistID, varID, varname);

          len = (int) strlen(varname);
702
703
          for (i = 0; i < len; i++)
            if (varname[i] == '-') break;
704

705
          if (i < len)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706
            for (j = i; j < len; j++) varname[j] = varname[j + 1];
707
708
709
710
711

          vlistInqVarLongname(vlistID, varID, varlongname);
          vlistInqVarUnits(vlistID, varID, varunits);
          fprintf(gdp, "%-15s", varname);

712
          if (nlev == 1) nlev = 0;
713
714
715

          fprintf(gdp, "  %3d", nlev);

716
          if (filetype == CDI_FILETYPE_GRB)
717
718
719
720
721
722
723
724
725
726
            {
              code = vlistInqVarCode(vlistID, varID);
              /*
              if      ( ltype == ZAXIS_SURFACE )  ltype = 1;
              else if ( ltype == ZAXIS_PRESSURE ) ltype = 99;
              else if ( nlev == 1 )  ltype = 1;
              else ltype = 99;
              */
              fprintf(gdp, "  %d,%d", code, ltype);
            }
727
          else if (filetype == CDI_FILETYPE_NC)
728
729
730
731
            {
              int xyz = vlistInqVarXYZ(vlistID, varID);

              fprintf(gdp, "  ");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
732
              if (vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT) fprintf(gdp, "t,");
733
              if (xyz == 321)
734
                {
735
                  if (nlev > 0) fprintf(gdp, "z,");
736
737
                  fprintf(gdp, "y,x");
                }
738
              else if (xyz == 312)
739
                {
740
                  if (nlev > 0) fprintf(gdp, "z,");
741
742
                  fprintf(gdp, "x,y");
                }
743
              else if (xyz == 231)
744
745
                {
                  fprintf(gdp, "y,");
746
                  if (nlev > 0) fprintf(gdp, "z,");
747
748
                  fprintf(gdp, "x");
                }
749
              else if (xyz == 132)
750
751
                {
                  fprintf(gdp, "x,");
752
                  if (nlev > 0) fprintf(gdp, "z,");
753
754
755
756
                  fprintf(gdp, "y");
                }
              else
                {
757
                  if (nlev > 0) fprintf(gdp, "z,");
758
759
760
761
762
763
                  fprintf(gdp, "y,x");
                }
            }
          else
            fprintf(gdp, "  99");

764
          if (varlongname[0] == 0)
765
766
767
768
            fprintf(gdp, "  %s", varname);
          else
            fprintf(gdp, "  %s", varlongname);

769
          if (varunits[0] != 0) fprintf(gdp, "  [%s]", varunits);
770
771
772

          fprintf(gdp, "\n");
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
773
774
775
776
777
    }

  fprintf(gdp, "ENDVARS\n");
}

778
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
779
write_map_grib1(const char *ctlfile, int map_version, int nrecords, int *intnum, float *fltnum, off_t *bignum)
780
781
782
{
  int i;
  struct gaindx indx;
783
  struct gaindxb indxb;
784
  FILE *mapfp;
785
  int hinum[5];
786

787
788
  memset(&indx, 0, sizeof(struct gaindx));

789
  mapfp = fopen(ctlfile, "w");
790
  if (mapfp == NULL) cdoAbort("Open failed on %s", ctlfile);
791

792
793
794
  indx.type = map_version;
  indx.hfnum = 0;
  if (map_version == 4)
795
    {
796
      indx.hinum = 5;
797
798
799
800
801
      indx.intnum = nrecords;
      indxb.bignum = 2 * nrecords;
    }
  else
    {
802
      indx.hinum = 4;
803
804
805
      indx.intnum = 3 * nrecords;
      indxb.bignum = 0;
    }
806
  indx.fltnum = 3 * nrecords;
807

808
809
  indx.hipnt = NULL;
  indx.hfpnt = NULL;
810
811
  indx.intpnt = NULL;
  indx.fltpnt = NULL;
812
  indxb.bigpnt = NULL;
813

814
  hinum[0] = map_version;
815
816
817
  hinum[1] = 1;
  hinum[2] = nrecords;
  hinum[3] = 255;
818
  hinum[4] = indxb.bignum;
819

820
  if (map_version == 2)
821
822
823
824
    {
      int nb, bcnt, rc, j;
      float fdum;
      unsigned char ibmfloat[4];
825

826
      /* calculate the size of the ver==1 index file */
827

Uwe Schulzweida's avatar
Uwe Schulzweida committed
828
      nb = 2 + (indx.hinum * 4) + /* version in byte 2, then 4 ints with number of each data type */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
           indx.hinum * sizeof(int) + indx.hfnum * sizeof(int) + indx.intnum * sizeof(int) + indx.fltnum * sizeof(float);
830

831
      /* add additional info */
832

833
      nb += 7;     /* base time (+ sec)  for compatibility with earlier version 2 maps */
834
      nb += 8 * 4; /* grvals for time <-> grid conversion */
835

836
      unsigned char *map = (unsigned char *) Malloc(nb);
837

838
839
      bcnt = 0;
      Put1Byte(map, bcnt, 0);
840
      Put1Byte(map, bcnt, map_version);
841

842
843
844
845
      Put4Byte(map, bcnt, indx.hinum);
      Put4Byte(map, bcnt, indx.hfnum);
      Put4Byte(map, bcnt, indx.intnum);
      Put4Byte(map, bcnt, indx.fltnum);
846

847
848
849
850
851
852
      Put2Byte(map, bcnt, 0); /* initial year   */
      Put1Byte(map, bcnt, 0); /* initial month  */
      Put1Byte(map, bcnt, 0); /* initial day    */
      Put1Byte(map, bcnt, 0); /* initial hour   */
      Put1Byte(map, bcnt, 0); /* initial minute */
      Put1Byte(map, bcnt, 0); /* initial second */
853

854
      if (indx.hinum)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
855
        for (i = 0; i < indx.hinum; i++) Put4Byte(map, bcnt, hinum[i]);
856

857
858
859
860
      if (indx.hfnum)
        {
          /* blank for now */
        }
861

862
      for (i = 0; i < nrecords; i++)
863
        {
864
865
          PutInt(map, bcnt, (int) bignum[i * 2]);
          PutInt(map, bcnt, (int) bignum[i * 2 + 1]);
866
867
          PutInt(map, bcnt, intnum[i]);
        }
868

869
      for (i = 0; i < indx.fltnum; i++)
870
        {
871
          fdum = fltnum[i];
872
          rc = flt2ibm(fdum, ibmfloat);
873
          if (rc < 0) cdoAbort("overflow in IBM float conversion");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
874
          for (j = 0; j < 4; j++) map[bcnt++] = ibmfloat[j];
875
876
877
878
        }

      /* write out the factors for converting from grid to absolute time */

879
      for (i = 0; i < 8; i++)
880
881
882
        {
          fdum = 0;
          rc = flt2ibm(fdum, ibmfloat);
883
          if (rc < 0) cdoAbort("overflow in IBM float conversion");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
884
          for (j = 0; j < 4; j++) map[bcnt++] = ibmfloat[j];
885
886
        }

887
      fwrite(map, 1, bcnt, mapfp);
888

889
      Free(map);
890
891
892
893
    }
  else
    {
      fwrite(&indx, sizeof(struct gaindx), 1, mapfp);
894
895
      if (indx.hinum > 0) fwrite(hinum, sizeof(int), indx.hinum, mapfp);
      if (map_version == 1)
896
        {
897
          std::vector<int> intnumbuf(indx.intnum);
898
          for (i = 0; i < nrecords; i++)
899
            {
900
901
902
              intnumbuf[i * 3 + 0] = (int) bignum[i * 2];
              intnumbuf[i * 3 + 1] = (int) bignum[i * 2 + 1];
              intnumbuf[i * 3 + 2] = intnum[i];
903
            }
904
          if (indx.intnum > 0) fwrite(intnumbuf.data(), sizeof(int), indx.intnum, mapfp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
905
          if (indx.fltnum > 0) fwrite(fltnum, sizeof(float), indx.fltnum, mapfp);
906
        }
907
      else
908
        {
909
          if (indx.intnum > 0) fwrite(intnum, sizeof(int), indx.intnum, mapfp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
910
911
          if (indx.fltnum > 0) fwrite(fltnum, sizeof(float), indx.fltnum, mapfp);
          if (indxb.bignum > 0) fwrite(bignum, sizeof(off_t), indxb.bignum, mapfp);
912
        }
913
    }
914

915
916
917
  fclose(mapfp);
}

918
/*
919
static
920
921
void write_map_grib2(const char *ctlfile, int map_version, int nrecords, int
*intnum, float *fltnum, off_t *bignum)
922
923
{
}
924
*/
925

926
927
void *
Gradsdes(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
928
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
929
  int gridID