stream_cgribex.c 63.4 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
#if HAVE_CONFIG_H
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
3
4
#endif

5
#include <limits.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
#include <stdio.h>

#include "dmemory.h"
#include "cdi.h"
10
#include "cdi_int.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
11
12
13
#include "file.h"
#include "varscan.h"
#include "datetime.h"
14
#include "stream_scan.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
15
#include "stream_grb.h"
16
#include "stream_cgribex.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
17

18
#ifdef HAVE_LIBCGRIBEX
19

Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
#include "cgribex.h"
21
22
23
24
25
26
27
28
29


typedef struct {
  int       sec0[2];
  int       sec1[1024];
  size_t    sec2len;
  int      *sec2;
  int       sec3[2];
  int       sec4[512];
30
31
  double   fsec2[512];
  double   fsec3[2];
32
33
34
}
cgribexrec_t;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36

typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
  int param;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
39
40
  int level1;
  int level2;
  int ltype;
41
  int tsteptype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
} compvar_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44


45
46
47
48
49
50
51
52
53
54
55
typedef struct
{
  void *gribbuffer;
  size_t gribbuffersize;
  unsigned char *pds;
  unsigned char *gds;
  unsigned char *bms;
  unsigned char *bds;
} cgribex_handle;


56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

static
void cgribexInit(cgribexrec_t *cgribexp)
{
  cgribexp->sec2len = 4096;
  cgribexp->sec2 = (int *) Malloc(cgribexp->sec2len*sizeof(int));
}


void *cgribexNew()
{
  cgribexrec_t *cgribexp = (cgribexrec_t *) Malloc(sizeof(cgribexrec_t));
  cgribexInit(cgribexp);
  return (void*)cgribexp;
}


void cgribexDelete(void *cgribex)
{
  cgribexrec_t *cgribexp = (cgribexrec_t *) cgribex;
  if (cgribexp)
    {
      if (cgribexp->sec2) Free(cgribexp->sec2);
      Free(cgribexp);
    }
}


84
85
86
int grib1Sections(unsigned char *gribbuffer, long gribbufsize, unsigned char **pdsp,
		  unsigned char **gdsp, unsigned char **bmsp, unsigned char **bdsp, long *gribrecsize);

87
size_t cgribexSection2Length(void *gribbuffer, size_t gribbuffersize)
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
{
  long sec2len = 0;

  if ( gribbuffersize && gribbuffer )
    {
      unsigned char *pds = NULL, *gds = NULL, *bms = NULL, *bds = NULL;
      long gribrecsize;
      int status = grib1Sections((unsigned char *)gribbuffer, (long)gribbuffersize, &pds, &gds, &bms, &bds, &gribrecsize);
      if (status >= 0 && gds != NULL) sec2len = (unsigned) ((gds[0] << 16) + (gds[1] << 8)  + (gds[3]));
    }

  return sec2len;
}


void *cgribex_handle_new_from_meassage(void *gribbuffer, size_t gribbuffersize)
{
  cgribex_handle *gh = (cgribex_handle*) Malloc(sizeof(cgribex_handle));
  gh->gribbuffer = NULL;
  gh->gribbuffersize = 0;
  gh->pds = NULL;

110
  if (gribbuffersize && gribbuffer)
111
112
113
114
    {
      unsigned char *pds = NULL, *gds = NULL, *bms = NULL, *bds = NULL;
      long gribrecsize;
      int status = grib1Sections((unsigned char *)gribbuffer, (long)gribbuffersize, &pds, &gds, &bms, &bds, &gribrecsize);
115
      if (status >= 0)
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
        {
          gh->gribbuffer = gribbuffer;
          gh->gribbuffersize = gribbuffersize;
          gh->pds = pds;
          gh->gds = gds;
          gh->bms = bms;
          gh->bds = bds;
        }
    }

  return (void*)gh;
}


void cgribex_handle_delete(void *gh)
{
132
  if (gh) Free(gh);
133
134
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
int cgribexGetGridType(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
{
138
  int gridtype = GRID_GENERIC;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139

140
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
  switch (ISEC2_GridType)
    {
143
    case  GRIB1_GTYPE_LATLON:     { gridtype = GRID_LONLAT;     break; }
144
    case  GRIB1_GTYPE_LATLON_ROT: { gridtype = GRID_PROJECTION; break; }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
    case  GRIB1_GTYPE_LCC:        { gridtype = CDI_PROJ_LCC;    break; }
146
    case  GRIB1_GTYPE_GAUSSIAN:   { gridtype = ISEC2_Reduced ? GRID_GAUSSIAN_REDUCED : GRID_GAUSSIAN; break; }
147
148
    case  GRIB1_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL;   break; }
    case  GRIB1_GTYPE_GME:        { gridtype = GRID_GME;        break; }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
    }
150
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151

Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
  return gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
155
static
156
bool cgribexGetIsRotated(int *isec2)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
{
158
  return (ISEC2_GridType == GRIB1_GTYPE_LATLON_ROT);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
static
162
bool cgribexGetZaxisHasBounds(int grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
{
164
  // clang-format off
165
  switch (grb_ltype)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
    {
167
    case GRIB1_LTYPE_SIGMA_LAYER:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
    case GRIB1_LTYPE_HYBRID_LAYER:
169
    case GRIB1_LTYPE_LANDDEPTH_LAYER: return true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
    }
171
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172

173
  return false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
int cgribexGetTimeUnit(int *isec1)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
{
179
  int timeunit = TUNIT_HOUR;
180
  static bool lprint = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
184
  switch ( ISEC1_TimeUnit )
    {
185
186
187
188
189
190
191
192
    case ISEC1_TABLE4_MINUTE:    timeunit = TUNIT_MINUTE;    break;
    case ISEC1_TABLE4_QUARTER:   timeunit = TUNIT_QUARTER;   break;
    case ISEC1_TABLE4_30MINUTES: timeunit = TUNIT_30MINUTES; break;
    case ISEC1_TABLE4_HOUR:      timeunit = TUNIT_HOUR;      break;
    case ISEC1_TABLE4_3HOURS:    timeunit = TUNIT_3HOURS;    break;
    case ISEC1_TABLE4_6HOURS:    timeunit = TUNIT_6HOURS;    break;
    case ISEC1_TABLE4_12HOURS:   timeunit = TUNIT_12HOURS;   break;
    case ISEC1_TABLE4_DAY:       timeunit = TUNIT_DAY;       break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
194
195
    default:
      if ( lprint )
	{
196
	  Message("GRIB time unit %d unsupported!", ISEC1_TimeUnit);
197
	  lprint = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
      break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
    }
201
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202

Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
  return timeunit;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
205
}

206
static
207
bool cgribexTimeIsFC(int *isec1)
208
{
209
  bool isFC = (ISEC1_TimeRange == 10 && ISEC1_TimePeriod1 == 0 && ISEC1_TimePeriod2 == 0) ? false : true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
  return isFC;
211
212
213
214
215
}

static
int cgribexGetTsteptype(int timerange)
{
216
  static bool lprint = true;
217

218
  // clang-format off
219
  int tsteptype = TSTEP_INSTANT;
220
221
222
223
224
225
226
227
228
229
230
231
  switch ( timerange )
    {
    case  0:  tsteptype = TSTEP_INSTANT;  break;
    case  1:  tsteptype = TSTEP_INSTANT2; break;
    case  2:  tsteptype = TSTEP_RANGE;    break;
    case  3:  tsteptype = TSTEP_AVG;      break;
    case  4:  tsteptype = TSTEP_ACCUM;    break;
    case  5:  tsteptype = TSTEP_DIFF;     break;
    case 10:  tsteptype = TSTEP_INSTANT3; break;
    default:
      if ( lprint )
	{
232
	  Message("Time range indicator %d unsupported, set to 0!", timerange);
233
	  lprint = false;
234
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
      break;
236
    }
237
  // clang-format on
238

Uwe Schulzweida's avatar
Uwe Schulzweida committed
239
  return tsteptype;
240
241
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, int iret)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
{
245
  bool compyinc = true;
246
  int gridtype = cgribexGetGridType(isec2);
247
  int projtype = (gridtype == GRID_PROJECTION && cgribexGetIsRotated(isec2)) ? CDI_PROJ_RLL : CDI_UNDEFID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
  if ( gridtype == CDI_PROJ_LCC )
249
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
      projtype = gridtype;
251
252
      gridtype = GRID_PROJECTION;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253

254
  if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED && iret != -801 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
    {
256
257
      int nlon = 0;
      for ( int ilat = 0; ilat < ISEC2_NumLat; ++ilat )
258
        if ( ISEC2_RowLon(ilat) > nlon ) nlon = ISEC2_RowLon(ilat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
      gridtype = GRID_GAUSSIAN;
260
261
      ISEC2_NumLon = nlon;
      ISEC4_NumValues = nlon*ISEC2_NumLat;
262
      compyinc = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263
264
    }

265
  grid_init(grid);
266
  cdiGridTypeInit(grid, gridtype, 0);
267

268
  if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_RLL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269
    {
270
271
      const bool ijDirectionIncrementGiven = gribbyte_get_bit(ISEC2_ResFlag, 1);
      const bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
272
      if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;
273

274
275
276
      const size_t nvalues = (size_t) ISEC4_NumValues;
      const size_t nlon = (size_t) ISEC2_NumLon;
      const size_t nlat = (size_t) ISEC2_NumLat;
277
278
279
280
281
282
283
      if ( nvalues != nlon*nlat )
        Error("numberOfPoints (%zu) and gridSize (%zu) differ!", nvalues, nlon*nlat);

      grid->size   = nvalues;
      grid->x.size = nlon;
      grid->y.size = nlat;

284
285
286
287
288
      if ( gridtype == GRID_GAUSSIAN ) grid->np = ISEC2_NumPar;
      grid->x.inc  = 0;
      grid->y.inc  = 0;
      grid->x.flag = 0;
      /* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
      {
290
291
292
        if ( grid->x.size > 1 )
          {
            bool recompinc = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293

294
295
296
297
298
299
300
            if ( ISEC2_LastLon < ISEC2_FirstLon )
              {
                if ( ISEC2_FirstLon >= 180000 )
                  ISEC2_FirstLon -= 360000;
                else
                  ISEC2_LastLon += 360000;
              }
301

302
            if ( ijDirectionIncrementGiven && ISEC2_LonIncr > 0 )
303
              {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
                if ( labs(ISEC2_LastLon - (ISEC2_FirstLon+ISEC2_LonIncr*((long)grid->x.size-1))) <= 2 )
305
                  {
306
307
                    recompinc = false;
                    grid->x.inc = ISEC2_LonIncr * 0.001;
308
                  }
309
              }
310

311
312
            /* recompute xinc if necessary */
            if ( recompinc ) grid->x.inc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid->x.size-1);
313

314
315
316
317
318
319
320
321
322
323
324
325
326
327
            /* correct xinc if necessary */
            if ( ISEC2_FirstLon == 0 && ISEC2_LastLon > 354000 && ISEC2_LastLon < 360000 )
              {
                double xinc = 360. / grid->x.size;
                if ( fabs(grid->x.inc-xinc) > 0.0 )
                  {
                    grid->x.inc = xinc;
                    if ( CDI_Debug ) Message("set xinc to %g", grid->x.inc);
                  }
              }
          }
        grid->x.first = ISEC2_FirstLon * 0.001;
        grid->x.last  = ISEC2_LastLon  * 0.001;
        grid->x.flag  = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
      }
329
330
      grid->y.flag = 0;
      /* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
      {
332
333
334
        if ( grid->y.size > 1 && compyinc )
          {
            bool recompinc = true;
335
            if ( ijDirectionIncrementGiven && ISEC2_LatIncr > 0 )
336
              {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
                if ( labs(ISEC2_LastLat - (ISEC2_FirstLat+ISEC2_LatIncr*((long)grid->y.size-1))) <= 2 )
338
339
340
341
342
343
344
345
346
347
348
349
                  {
                    recompinc = false;
                    grid->y.inc = ISEC2_LatIncr * 0.001;
                  }
              }

            /* recompute yinc if necessary */
            if ( recompinc ) grid->y.inc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid->y.size - 1);
          }
        grid->y.first = ISEC2_FirstLat * 0.001;
        grid->y.last  = ISEC2_LastLat  * 0.001;
        grid->y.flag  = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
      }
351
352
353
    }
  else if ( gridtype == GRID_GAUSSIAN_REDUCED )
    {
354
355
      const bool ijDirectionIncrementGiven = gribbyte_get_bit(ISEC2_ResFlag, 1);
      const bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
356
      if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;
357
      grid->np      = ISEC2_NumPar;
358
      grid->size    = (size_t)ISEC4_NumValues;
359
      grid->rowlon  = ISEC2_RowLonPtr;
360
361
      grid->nrowlon = (size_t)ISEC2_NumLat;
      grid->y.size  = (size_t)ISEC2_NumLat;
362
363
364
365
      grid->x.inc   = 0;
      grid->y.inc   = 0;
      grid->x.flag  = 0;
      /* if ( ISEC2_FirstLon != 0 || ISEC2_LastLon != 0 ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
      {
367
368
369
370
        if ( grid->x.size > 1 )
          {
            if ( ISEC2_LastLon < ISEC2_FirstLon && ISEC2_LastLon < 0 ) ISEC2_LastLon += 360000;

371
            if ( ijDirectionIncrementGiven && ISEC2_LonIncr > 0 )
372
373
374
375
376
377
378
              grid->x.inc = ISEC2_LonIncr * 0.001;
            else
              grid->x.inc = (ISEC2_LastLon - ISEC2_FirstLon) * 0.001 / (grid->x.size - 1);
          }
        grid->x.first = ISEC2_FirstLon * 0.001;
        grid->x.last  = ISEC2_LastLon  * 0.001;
        grid->x.flag  = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
      }
380
381
      grid->y.flag = 0;
      /* if ( ISEC2_FirstLat != 0 || ISEC2_LastLat != 0 ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
382
      {
383
384
        if ( grid->y.size > 1 )
          {
385
            if ( ijDirectionIncrementGiven && ISEC2_LatIncr > 0 )
386
387
388
389
390
391
392
              grid->y.inc = ISEC2_LatIncr * 0.001;
            else
              grid->y.inc = (ISEC2_LastLat - ISEC2_FirstLat) * 0.001 / (grid->y.size - 1);
          }
        grid->y.first = ISEC2_FirstLat * 0.001;
        grid->y.last  = ISEC2_LastLat  * 0.001;
        grid->y.flag  = 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393
394
      }
    }
395
  else if ( projtype == CDI_PROJ_LCC )
396
    {
397
      const bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
398
399
      if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;

400
401
402
      const size_t nvalues = (size_t) ISEC4_NumValues;
      const size_t nlon = (size_t) ISEC2_NumLon;
      const size_t nlat = (size_t) ISEC2_NumLat;
403
404
      if ( nvalues != nlon*nlat )
        Error("numberOfPoints (%zu) and gridSize (%zu) differ!", nvalues, nlon*nlat);
405

406
407
408
      grid->size   = nvalues;
      grid->x.size = nlon;
      grid->y.size = nlat;
409

410
411
412
413
414
415
416
417
      grid->x.first = 0;
      grid->x.last  = 0;
      grid->x.inc   = ISEC2_Lambert_dx;
      grid->y.first = 0;
      grid->y.last  = 0;
      grid->y.inc   = ISEC2_Lambert_dy;
      grid->x.flag  = 2;
      grid->y.flag  = 2;
418
419
420
    }
  else if ( gridtype == GRID_SPECTRAL )
    {
421
      grid->size  = (size_t) ISEC4_NumValues;
422
      grid->trunc = ISEC2_PentaJ;
423
424
      grid->lcomplex = (ISEC2_RepMode == 2) ? 1 : 0;
   }
425
426
  else if ( gridtype == GRID_GME )
    {
427
      grid->size  = (size_t) ISEC4_NumValues;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
429
430
431
      grid->gme.nd  = ISEC2_GME_ND;
      grid->gme.ni  = ISEC2_GME_NI;
      grid->gme.ni2 = ISEC2_GME_NI2;
      grid->gme.ni3 = ISEC2_GME_NI3;
432
433
434
    }
  else if ( gridtype == GRID_GENERIC )
    {
435
      grid->size  = (size_t) ISEC4_NumValues;
436
437
438
439
440
441
442
      grid->x.size = 0;
      grid->y.size = 0;
    }
  else
    {
      Error("Unsupported grid type: %s", gridNamePtr(gridtype));
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443

444
445
  grid->type = gridtype;
  grid->projtype = projtype;
446
447
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
449
450
451
452
453
454
455
456
457
static
void cgribexGetLevel(int *isec1, int *leveltype, int *level1, int *level2)
{
  *leveltype = ISEC1_LevelType;
  *level1 = ISEC1_Level1;
  *level2 = ISEC1_Level2;
  if ( *leveltype == GRIB1_LTYPE_ISOBARIC ) *level1 *= 100;
  else if ( *leveltype == GRIB1_LTYPE_99 || *leveltype == GRIB1_LTYPE_ISOBARIC_PA ) *leveltype = GRIB1_LTYPE_ISOBARIC;
}

458
static
459
void cgribexAddRecord(stream_t *streamptr, cgribexrec_t *cgribexp, int param, size_t recsize,
460
                      off_t position, int comptype, int lmv, int iret)
461
{
462
  int *isec1 = cgribexp->sec1;
463
  int *isec2 = cgribexp->sec2;
464
  int *isec4 = cgribexp->sec4;
465
466
467
  double *fsec2 = cgribexp->fsec2;
  double *fsec3 = cgribexp->fsec3;

468
469
  int datatype = (ISEC4_NumBits > 0 && ISEC4_NumBits <= 32) ? ISEC4_NumBits : CDI_DATATYPE_PACK;

470
  int varID;
471
472
  int levelID = 0;

473
474
475
  const int vlistID = streamptr->vlistID;
  const int tsID    = streamptr->curTsID;
  const int recID   = recordNewEntry(streamptr, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476
  record_t *record = &streamptr->tsteps[tsID].records[recID];
477

478
479
  const int tsteptype = cgribexGetTsteptype(ISEC1_TimeRange);
  const int numavg    = ISEC1_AvgNum;
480

Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
482
  int leveltype, level1, level2;
  cgribexGetLevel(isec1, &leveltype, &level1, &level2);
483

Uwe Schulzweida's avatar
Uwe Schulzweida committed
484
  /* fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, leveltype); */
485

486
  record->size      = recsize;
487
488
489
490
  record->position  = position;
  record->param     = param;
  record->ilevel    = level1;
  record->ilevel2   = level2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
491
  record->ltype     = leveltype;
Thomas Jahns's avatar
Thomas Jahns committed
492
  record->tsteptype = (short)tsteptype;
493

Uwe Schulzweida's avatar
Uwe Schulzweida committed
494
  grid_t *gridptr = (grid_t*) Malloc(sizeof(*gridptr));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
  cgribexGetGrid(streamptr, isec2, isec4, gridptr, iret);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
496

Uwe Schulzweida's avatar
Uwe Schulzweida committed
497
  struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, gridptr, 0);
498
  int gridID = gridAdded.Id;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
499
500
501
502
  if ( gridAdded.isNew )
    {
      if ( gridptr->nrowlon )
        {
503
          const size_t nrowlon = (size_t) gridptr->nrowlon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
505
506
507
          int *rowlon = gridptr->rowlon;
          gridptr->rowlon = (int*) Malloc(nrowlon * sizeof(int));
          memcpy(gridptr->rowlon, rowlon, nrowlon * sizeof(int));
        }
508
509
      else if ( gridptr->projtype == CDI_PROJ_RLL )
        {
510
511
512
          const double xpole =   ISEC2_LonSP*0.001 - 180;
          const double ypole = - ISEC2_LatSP*0.001;
          const double angle = - FSEC2_RotAngle;
513
          gridDefParamRLL(gridID, xpole, ypole, angle);
514
        }
515
516
      else if ( gridptr->projtype == CDI_PROJ_LCC )
        {
517
518
519
          double a = 6367470., rf = 0;
          bool earthIsOblate = gribbyte_get_bit(ISEC2_ResFlag, 2);
          if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
520
521
522
          const double xval_0 = ISEC2_FirstLon * 0.001;
          const double yval_0 = ISEC2_FirstLat * 0.001;
          const double lon_0  = ISEC2_Lambert_Lov * 0.001;
523
524
          double lat_1  = ISEC2_Lambert_LatS1 * 0.001;
          double lat_2  = ISEC2_Lambert_LatS2 * 0.001;
525
          const bool lsouth = gribbyte_get_bit(ISEC2_Lambert_ProjFlag, 1);
526
          if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
527

528
          const double lat_0 = lat_2;
529
530
          double x_0 = CDI_grid_missval;
          double y_0 = CDI_grid_missval;
531
532
533
534

          if ( proj_lonlat_to_lcc_func )
            {
              x_0 = xval_0; y_0 = yval_0;
535
536
              proj_lonlat_to_lcc_func(CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
              if ( IS_NOT_EQUAL(x_0, CDI_grid_missval) && IS_NOT_EQUAL(y_0, CDI_grid_missval) )
537
538
                { x_0 = -x_0; y_0 = -y_0; }
            }
539
          gridDefParamLCC(gridID, CDI_grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0);
540
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
541
542
543
    }
  else
    Free(gridptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
544

545
  const int zaxistype = grib1ltypeToZaxisType(leveltype);
546
  if ( zaxistype == ZAXIS_HYBRID || zaxistype == ZAXIS_HYBRID_HALF )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547
    {
548
      const size_t vctsize = (size_t)ISEC2_NumVCP;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
551
552
      double *vctptr = &fsec2[10];
      varDefVCT(vctsize, vctptr);
    }

553
  const bool lbounds = cgribexGetZaxisHasBounds(leveltype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
554

555
556
  if ( datatype > 32 ) datatype = CDI_DATATYPE_PACK32;
  if ( datatype <  0 ) datatype = CDI_DATATYPE_PACK;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557

558
  varAddRecord(recID, param, gridID, zaxistype, lbounds, level1, level2, 0, 0,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
559
	       datatype, &varID, &levelID, tsteptype, numavg, leveltype, -1,
560
               NULL, NULL, NULL, NULL, NULL, NULL, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
561

Thomas Jahns's avatar
Thomas Jahns committed
562
563
  record->varID   = (short)varID;
  record->levelID = (short)levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
564

565
  varDefCompType(varID, comptype);
566

567
568
569
  if ( ISEC1_LocalFLag )
    {
      if      ( ISEC1_CenterID == 78  && isec1[36] == 253 ) // DWD local extension
570
571
572
573
574
        {
          varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, isec1[52]);
          varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, isec1[53]);
          varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, isec1[54]);
        }
575
      else if ( ISEC1_CenterID == 252 && isec1[36] ==   1 ) // MPIM local extension
576
577
578
579
580
        {
          varDefKeyInt(varID, CDI_KEY_TYPEOFENSEMBLEFORECAST, isec1[37]);
          varDefKeyInt(varID, CDI_KEY_NUMBEROFFORECASTSINENSEMBLE, isec1[39]);
          varDefKeyInt(varID, CDI_KEY_PERTURBATIONNUMBER, isec1[38]);
        }
581
    }
582

583
584
  if ( lmv ) varDefMissval(varID, FSEC3_MissVal);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
585
  if ( varInqInst(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
    {
587
588
      const int center    = ISEC1_CenterID;
      const int subcenter = ISEC1_SubCenterID;
589
      int instID    = institutInq(center, subcenter, NULL, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
590
      if ( instID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
591
592
593
594
	instID = institutDef(center, subcenter, NULL, NULL);
      varDefInst(varID, instID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
595
  if ( varInqModel(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
596
    {
597
      int modelID = modelInq(varInqInst(varID), ISEC1_ModelID, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
      if ( modelID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
599
600
601
602
	modelID = modelDef(varInqInst(varID), ISEC1_ModelID, NULL);
      varDefModel(varID, modelID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
603
  if ( varInqTable(varID) == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604
    {
605
      int tableID = tableInq(varInqModel(varID), ISEC1_CodeTable, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
606
      if ( tableID == CDI_UNDEFID )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
608
609
610
	tableID = tableDef(varInqModel(varID), ISEC1_CodeTable, NULL);
      varDefTable(varID, tableID);
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
612
  streamptr->tsteps[tsID].nallrecs++;
  streamptr->nrecs++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613
614
}

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
static
void MCH_get_undef(int *isec1, double *undef_pds, double *undef_eps)
{
  /* 2010-01-13: Oliver Fuhrer */
  if ( ISEC1_CenterID == 215 ) {
    if (isec1[34] != 0 && isec1[34] != 255) {
      if (isec1[34] & 2) {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,-isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,-isec1[35]);
        }
        *undef_eps = pow(10.0,-isec1[35]-1);
      } else {
        if (isec1[34] & 1) {
          *undef_pds = -0.99*pow(10.0,+isec1[35]);
        } else {
          *undef_pds = +0.99*pow(10.0,+isec1[35]);
        }
        *undef_eps = pow(10.0,isec1[35]-1);
      }
    }
  }
}

static
641
void cgribexDecodeHeader(cgribexrec_t *cgribexp, int *gribbuffer, int recsize, int *lmv, int *iret)
642
{
643
  int *isec0 = cgribexp->sec0;
644
  int *isec1 = cgribexp->sec1;
645
  int *isec2 = cgribexp->sec2;
646
  int *isec3 = cgribexp->sec3;
647
  int *isec4 = cgribexp->sec4;
648
649
650
  double *fsec2 = cgribexp->fsec2;
  double *fsec3 = cgribexp->fsec3;

651
  int ipunp = 0, iword = 0;
652

653
  memset(isec1, 0, 256*sizeof(int));
654
  memset(isec2, 0, 32*sizeof(int));
655

656
  double *fsec4 = NULL;
657
  gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
658
	   ipunp, (int *) gribbuffer, recsize, &iword, "J", iret);
659

Uwe Schulzweida's avatar
Uwe Schulzweida committed
660
661
  if ( !(ISEC1_Sec2Or3Flag & 128) ) isec2[0] = -1; // default generic grid

662
663
664
665
666
667
668
669
670
671
  *lmv = 0;

  if ( ISEC1_CenterID == 215 && (isec1[34] != 0 && isec1[34] != 255) )
    {
      double undef_pds, undef_eps;
      MCH_get_undef(isec1, &undef_pds, &undef_eps);
      FSEC3_MissVal = undef_pds;
      *lmv = 1;
    }
}
672
673

static
674
compvar_t cgribexVarSet(int param, int level1, int level2, int leveltype, int trange)
675
{
676
  int tsteptype = cgribexGetTsteptype(trange);
677

678
  compvar_t compVar;
679
680
681
682
683
  compVar.param     = param;
  compVar.level1    = level1;
  compVar.level2    = level2;
  compVar.ltype     = leveltype;
  compVar.tsteptype = tsteptype;
684

Uwe Schulzweida's avatar
Uwe Schulzweida committed
685
  return compVar;
686
687
}

688
689
static inline
int cgribexVarCompare(compvar_t compVar, record_t record, int flag)
690
{
691
692
693
694
  bool vinst = compVar.tsteptype == TSTEP_INSTANT || compVar.tsteptype == TSTEP_INSTANT2 || compVar.tsteptype == TSTEP_INSTANT3;
  bool rinst = record.tsteptype == TSTEP_INSTANT || record.tsteptype == TSTEP_INSTANT2 || record.tsteptype == TSTEP_INSTANT3;
  int tstepDiff = (!((flag == 0) & (vinst && rinst)))
                & (compVar.tsteptype != record.tsteptype);
Thomas Jahns's avatar
Thomas Jahns committed
695
696
697
698
699
  int rstatus = (compVar.param != record.param)
    |           (compVar.level1 != record.ilevel)
    |           (compVar.level2 != record.ilevel2)
    |           (compVar.ltype != record.ltype)
    |           tstepDiff;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
700
  return rstatus;
701
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
702

703
704
705
#define gribWarning(text, nrecs, timestep, paramstr, level1, level2) \
            Warning("Record %2d (id=%s lev1=%d lev2=%d) timestep %d: %s", nrecs, paramstr, level1, level2, timestep, text)

706
707
static inline
void cgribexScanTsFixNtsteps(stream_t *streamptr, off_t recpos)
708
709
710
{
  if ( streamptr->ntsteps == -1 )
    {
711
      const int tsID = tstepsNewEntry(streamptr);
712
713
714
      if ( tsID != streamptr->rtsteps )
	Error("Internal error. tsID = %d", tsID);

715
      streamptr->tsteps[tsID-1].next   = true;
716
717
718
719
      streamptr->tsteps[tsID].position = recpos;
    }
}

720
721
722
723
724
725
726
727
728
729
730
731
732
static
void cgribexSkipRecords(const int fileID)
{
  int nskip = cdiSkipRecords;
  while (nskip-- > 0)
    {
      const size_t recsize = gribGetSize(fileID);
      if (recsize == 0) Error("Skipping of %d records failed!", cdiSkipRecords);

      const off_t recpos = fileGetPos(fileID);
      fileSetPos(fileID, recpos, SEEK_CUR);
    }
}
733

Uwe Schulzweida's avatar
Uwe Schulzweida committed
734
int cgribexScanTimestep1(stream_t *streamptr)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
735
{
736
  int lmv = 0, iret = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
737
  off_t recpos = 0;
738
  void *gribbuffer = NULL;
739
  size_t buffersize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
740
  int leveltype = 0, level1 = 0, level2 = 0, vdate = 0, vtime = 0;
741
  DateTime datetime, datetime0 = { LONG_MIN, LONG_MIN };
742
  unsigned recID;
743
  int nrecs_scanned = 0;
744
745
  bool warn_time = true;
  bool warn_numavg = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
746
  int taxisID = -1;
747
748
  int rdate = 0, rtime = 0, tunit = 0;
  bool fcast = false;
749
  char paramstr[32];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
750
751

  streamptr->curTsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
752

753
  cgribexrec_t *cgribexp = (cgribexrec_t *) streamptr->record->cgribexp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
754

755
  const int tsID = tstepsNewEntry(streamptr);
756
  if (tsID != 0) Error("Internal problem! tstepsNewEntry returns %d", tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
757

758
  taxis_t *taxis = &streamptr->tsteps[tsID].taxis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759

760
  const int fileID = streamptr->fileID;
761

762
  if (cdiSkipRecords) cgribexSkipRecords(fileID);
763

764
  unsigned nrecs = 0;
765
  while (true)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
766
    {
767
      size_t recsize = gribGetSize(fileID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
768
769
770
771
      recpos  = fileGetPos(fileID);

      if ( recsize == 0 )
	{
772
	  if ( nrecs == 0 ) Error("No GRIB records found!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
773
	  streamptr->ntsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
774
775
	  break;
	}
776

777
      ensureBufferSize(recsize, &buffersize, &gribbuffer);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
778

779
780
781
      size_t readsize = recsize;
      // Search for next 'GRIB', read the following record, and position file offset after it.
      if (gribRead(fileID, gribbuffer, &readsize)) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
782

783
      const int comptype = grbDecompress(recsize, &buffersize, &gribbuffer);
784

785
      const size_t sec2len = cgribexSection2Length(gribbuffer, buffersize);
786
787
788
789
790
791
792
      if (sec2len > cgribexp->sec2len)
        {
          cgribexp->sec2len = sec2len;
          cgribexp->sec2 = (int *) Realloc(cgribexp->sec2, sec2len*sizeof(int));
        }

      int *isec1 = cgribexp->sec1;
793

794
      nrecs_scanned++;
795
      cgribexDecodeHeader(cgribexp, (int *) gribbuffer, (int)recsize, &lmv, &iret);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
796

797
      const int param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
798
799
      cdiParamToString(param, paramstr, sizeof(paramstr));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
800
      cgribexGetLevel(isec1, &leveltype, &level1, &level2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
801

Uwe Schulzweida's avatar
Uwe Schulzweida committed
802
803
      gribDateTime(isec1, &vdate, &vtime);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
804
805
806
807
808
809
      if ( nrecs == 0 )
	{
	  datetime0.date = vdate;
	  datetime0.time = vtime;
	  rdate = gribRefDate(isec1);
	  rtime = gribRefTime(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
810
	  tunit = cgribexGetTimeUnit(isec1);
811
	  fcast = cgribexTimeIsFC(isec1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
812
813
814
	}
      else
	{
815
816
	  datetime.date = vdate;
	  datetime.time = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
817
	  compvar_t compVar = cgribexVarSet(param, level1, level2, leveltype, ISEC1_TimeRange);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
818
819
	  for ( recID = 0; recID < nrecs; recID++ )
	    {
820
	      if ( cgribexVarCompare(compVar, streamptr->tsteps[0].records[recID], 0) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
821
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
822

823
	  if ( CDI_inventory_mode == 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
824
825
826
	    {
	      if ( recID < nrecs ) break;
	      if ( warn_time )
827
		if ( datetimeCmp(datetime, datetime0) != 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
828
		  {
829
                    gribWarning("Inconsistent verification time!", nrecs_scanned, tsID+1, paramstr, level1, level2);
830
		    warn_time = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
831
832
833
834
		  }
	    }
	  else
	    {
835
	      if ( datetimeCmp(datetime, datetime0) != 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
836
837
838

	      if ( recID < nrecs )
		{
839
		  gribWarning("Parameter already exist, skipped!", nrecs_scanned, tsID+1, paramstr, level1, level2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
840
841
842
		  continue;
		}
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
843
844
845
846
847
848
	}

      if ( ISEC1_AvgNum )
	{
	  if (  taxis->numavg && warn_numavg && (taxis->numavg != ISEC1_AvgNum) )
	    {
849
	      Warning("Changing numavg from %d to %d not supported!", taxis->numavg, ISEC1_AvgNum);
850
	      warn_numavg = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
851
852
853
854
855
856
857
858
859
860
	    }
	  else
	    {
	      taxis->numavg = ISEC1_AvgNum;
	    }
	}

      nrecs++;

      if ( CDI_Debug )
861
	Message("Read record %2d (id=%s lev1=%d lev2=%d) %8d %6d", nrecs_scanned, paramstr, level1, level2, vdate, vtime);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
862

863
      cgribexAddRecord(streamptr, cgribexp, param, recsize, recpos, comptype, lmv, iret);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
864
865
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
866
  streamptr->rtsteps = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
867

Uwe Schulzweida's avatar
Uwe Schulzweida committed
868
  if ( nrecs == 0 ) return CDI_EUFSTRUCT;
869

870
  cdi_generate_vars(streamptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
871

872
873
874
875
876
  taxis->unit = tunit;
  taxis->type = fcast ? TAXIS_RELATIVE : TAXIS_ABSOLUTE;
  taxisID = taxisCreate(taxis->type);
  taxis->rdate = rdate;
  taxis->rtime = rtime;
877
878
  taxis->vdate = (int)datetime0.date;
  taxis->vtime = (int)datetime0.time;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
879

880
  int vlistID = streamptr->vlistID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
881
882
  vlistDefTaxis(vlistID, taxisID);

883
  unsigned nrecords = (unsigned)streamptr->tsteps[0].nallrecs;
884
  if ( nrecords < (unsigned)streamptr->tsteps[0].recordSize )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
885
    {
886
      streamptr->tsteps[0].recordSize = (int)nrecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
887
      streamptr->tsteps[0].records =
888
        (record_t *) Realloc(streamptr->tsteps[0].records, nrecords*sizeof(record_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
889
890
    }

891
  streamptr->tsteps[0].recIDs = (int *) Malloc(nrecords*sizeof(int));
892
  streamptr->tsteps[0].nrecs = (int)nrecords;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
893
  for ( recID = 0; recID < nrecords; recID++ )
894
    streamptr->tsteps[0].recIDs[recID] = (int)recID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
895

Uwe Schulzweida's avatar
Uwe Schulzweida committed
896
  streamptr->record->buffer     = gribbuffer;
897
  streamptr->record->buffersize = buffersize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
898

899
  cgribexScanTsFixNtsteps(streamptr, recpos);
900
  streamScanTimeConstAdjust(streamptr, taxis);
901

Uwe Schulzweida's avatar
Uwe Schulzweida committed
902
  return</