cdf_read.c 31.4 KB
Newer Older
1 2
#ifdef HAVE_CONFIG_H
#include "config.h"
3 4 5 6 7 8 9 10 11 12
#endif

#ifdef HAVE_LIBNETCDF

#include <limits.h>
#include <float.h>

#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
13
#include "stream_grb.h"
14 15 16
#include "stream_cdf.h"
#include "cdf_int.h"
#include "vlist.h"
17
#include "vlist_var.h"
18 19 20 21 22


static
void cdfReadGridTraj(stream_t *streamptr, int gridID)
{
23 24
  const int vlistID = streamptr->vlistID;
  const int fileID  = streamptr->fileID;
25

26 27 28
  const int gridindex = vlistGridIndex(vlistID, gridID);
  const int lonID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_X];
  const int latID = streamptr->ncgrid[gridindex].ncIDs[CDF_VARID_Y];
29

30 31
  const int tsID = streamptr->curTsID;
  const size_t index = (size_t)tsID;
32 33 34 35 36 37 38 39 40 41 42 43

  double xlon, xlat;
  cdf_get_var1_double(fileID, lonID, &index, &xlon);
  cdf_get_var1_double(fileID, latID, &index, &xlat);

  gridDefXvals(gridID, &xlon);
  gridDefYvals(gridID, &xlat);
}

static
void cdfGetSlapDescription(stream_t *streamptr, int varID, size_t (*start)[4], size_t (*count)[4])
{
44 45 46 47 48 49
  const int vlistID = streamptr->vlistID;
  const int tsID = streamptr->curTsID;
  const int gridID = vlistInqVarGrid(vlistID, varID);
  const int zaxisID = vlistInqVarZaxis(vlistID, varID);
  const int timetype = vlistInqVarTimetype(vlistID, varID);
  const int gridindex = vlistGridIndex(vlistID, gridID);
50 51 52

  if ( CDI_Debug ) Message("tsID = %d", tsID);

53
  int xid = CDI_UNDEFID, yid = CDI_UNDEFID;
54 55 56 57 58 59
  if ( gridInqType(gridID) == GRID_TRAJECTORY )
    {
      cdfReadGridTraj(streamptr, gridID);
    }
  else
    {
60
      xid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
61
      yid = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_Y];
62
    }
63 64
  const int zaxisindex = vlistZaxisIndex(vlistID, zaxisID);
  const int zid = streamptr->zaxisID[zaxisindex];
65 66 67 68 69 70 71 72

  int ndims = 0;
#define addDimension(startCoord, length) do \
    { \
      (*start)[ndims] = startCoord; \
      (*count)[ndims] = length; \
      ndims++; \
    } while(0)
73
  if ( timetype != TIME_CONSTANT ) addDimension((size_t)tsID, 1);
74
  if ( zid != CDI_UNDEFID ) addDimension(0, (size_t)zaxisInqSize(zaxisID));
75 76
  if ( yid != CDI_UNDEFID ) addDimension(0, gridInqYsize(gridID));
  if ( xid != CDI_UNDEFID ) addDimension(0, gridInqXsize(gridID));
77 78 79 80 81 82 83 84 85 86
#undef addDimension

  assert(ndims <= (int)(sizeof(*start)/sizeof(**start)));
  assert(ndims <= (int)(sizeof(*count)/sizeof(**count)));

  if ( CDI_Debug )
    for (int idim = 0; idim < ndims; idim++)
      Message("dim = %d  start = %d  count = %d", idim, start[idim], count[idim]);
}

87 88
// Scans the data array for missVals, optionally applying first a scale factor and then an offset.
// Returns the number of missing + out-of-range values encountered.
89
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
size_t cdfDoInputDataTransformationDP(int vlistID, int varID, size_t valueCount, double *data)
91
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
  double missVal = vlistInqVarMissval(vlistID, varID);
93
  const int haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94 95 96 97 98 99 100 101
  double validRange[2];
  if (!(haveMissVal && vlistInqVarValidrange(vlistID, varID, validRange)))
    validRange[0] = DBL_MIN, validRange[1] = DBL_MAX;
  double validMin = validRange[0];
  double validMax = validRange[1];
  double offset = vlistInqVarAddoffset(vlistID, varID);
  double scaleFactor = vlistInqVarScalefactor(vlistID, varID);

102 103 104
  const bool missValIsNaN = DBL_IS_NAN(missVal);
  const int haveOffset = IS_NOT_EQUAL(offset, 0.0);
  const int haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1.0);
105 106 107 108 109
  size_t missValCount = 0;

  if ( IS_EQUAL(validMin, VALIDMISS) ) validMin = DBL_MIN;
  if ( IS_EQUAL(validMax, VALIDMISS) ) validMax = DBL_MAX;

110
  const int haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin, DBL_MIN));
111 112
  assert(!haveRangeCheck || haveMissVal);

113
  switch (haveMissVal | (haveScaleFactor << 1) | (haveOffset << 2) | (haveRangeCheck << 3))
114
    {
115
    case 15: // haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset
116 117 118 119 120 121 122 123 124
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? missVal
            : isMissVal ? data[i] : data[i] * scaleFactor + offset;
        }
      break;
125
    case 13: // haveRangeCheck & haveMissVal & haveOffset
126 127 128 129 130 131 132 133 134
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? missVal
            : isMissVal ? data[i] : data[i] + offset;
        }
      break;
135
    case 11: // haveRangeCheck & haveMissVal & haveScaleFactor
136 137 138 139 140 141 142 143 144
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? missVal
            : isMissVal ? data[i] : data[i] * scaleFactor;
        }
      break;
145
    case 9: // haveRangeCheck & haveMissVal
146 147 148 149 150 151 152 153
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? missVal : data[i];
        }
      break;
154
    case 7: // haveMissVal & haveScaleFactor & haveOffset
155 156 157 158 159 160
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] = data[i] * scaleFactor + offset;
      break;
161
    case 6: // haveOffset & haveScaleFactor
162 163 164
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] = data[i] * scaleFactor + offset;
      break;
165
    case 5: // haveMissVal & haveOffset
166 167 168 169 170 171
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] += offset;
      break;
172
    case 4: // haveOffset
173 174 175
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] += offset;
      break;
176
    case 3: // haveMissVal & haveScaleFactor
177 178 179 180 181 182
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] *= scaleFactor;
      break;
183
    case 2: // haveScaleFactor
184 185 186
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] *= scaleFactor;
      break;
187 188 189 190 191 192 193 194 195 196 197
    case 1: // haveMissVal
      if (missValIsNaN)
        {
          for ( size_t i = 0; i < valueCount; i++ )
            missValCount += (size_t)DBL_IS_NAN(data[i]);
        }
      else
        {
          for ( size_t i = 0; i < valueCount; i++ )
            missValCount += (size_t)IS_EQUAL(data[i], missVal);
        }
198 199 200 201 202 203 204
      break;
    }

  return missValCount;
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
size_t cdfDoInputDataTransformationSP(int vlistID, int varID, size_t valueCount, float *data)
206
 {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
  double missVal = vlistInqVarMissval(vlistID, varID);
208
  const int haveMissVal = vlistInqVarMissvalUsed(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209 210 211 212 213 214 215 216
  double validRange[2];
  if (!(haveMissVal && vlistInqVarValidrange(vlistID, varID, validRange)))
    validRange[0] = DBL_MIN, validRange[1] = DBL_MAX;
  double validMin = validRange[0];
  double validMax = validRange[1];
  double offset = vlistInqVarAddoffset(vlistID, varID);
  double scaleFactor = vlistInqVarScalefactor(vlistID, varID);

217 218 219
  const bool missValIsNaN = DBL_IS_NAN(missVal);
  const int haveOffset = IS_NOT_EQUAL(offset, 0.0);
  const int haveScaleFactor = IS_NOT_EQUAL(scaleFactor, 1.0);
220 221 222 223 224
  size_t missValCount = 0;

  if ( IS_EQUAL(validMin, VALIDMISS) ) validMin = DBL_MIN;
  if ( IS_EQUAL(validMax, VALIDMISS) ) validMax = DBL_MAX;

225
  const int haveRangeCheck = (IS_NOT_EQUAL(validMax, DBL_MAX)) | (IS_NOT_EQUAL(validMin, DBL_MIN));
226 227
  assert(!haveRangeCheck || haveMissVal);

228
  switch (haveMissVal | (haveScaleFactor << 1) | (haveOffset << 2) | (haveRangeCheck << 3))
229
    {
230
    case 15: // haveRangeCheck & haveMissVal & haveScaleFactor & haveOffset
231 232 233 234 235 236 237 238 239
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? (float)missVal
            : isMissVal ? data[i] : (float)(data[i] * scaleFactor + offset);
        }
      break;
240
    case 13: // haveRangeCheck & haveMissVal & haveOffset
241 242 243 244 245 246 247 248 249
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? (float)missVal
            : isMissVal ? data[i] : (float)(data[i] + offset);
        }
      break;
250
    case 11: // haveRangeCheck & haveMissVal & haveScaleFactor
251 252 253 254 255 256 257 258 259
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? (float)missVal
            : isMissVal ? data[i] : (float)(data[i] * scaleFactor);
        }
      break;
260
    case 9: // haveRangeCheck & haveMissVal
261 262 263 264 265 266 267 268
      for ( size_t i = 0; i < valueCount; i++ )
        {
          int outOfRange = data[i] < validMin || data[i] > validMax;
          int isMissVal = DBL_IS_EQUAL(data[i], missVal);
          missValCount += (size_t)(outOfRange | isMissVal);
          data[i] = outOfRange ? (float)missVal : data[i];
        }
      break;
269
    case 7: // haveMissVal & haveScaleFactor & haveOffset
270 271 272 273 274 275
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] = (float)(data[i] * scaleFactor + offset);
      break;
276
    case 6: // haveOffset & haveScaleFactor
277 278 279
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] = (float)(data[i] * scaleFactor + offset);
      break;
280
    case 5: // haveMissVal & haveOffset
281 282 283 284 285 286
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] = (float)(data[i] + offset);
      break;
287
    case 4: // haveOffset
288 289 290
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] = (float)(data[i] + offset);
      break;
291
    case 3: // haveMissVal & haveScaleFactor
292 293 294 295 296 297
      for ( size_t i = 0; i < valueCount; i++ )
        if ( DBL_IS_EQUAL(data[i], missVal) )
          missValCount++;
        else
          data[i] = (float)(data[i] * scaleFactor);
      break;
298
    case 2: // haveScaleFactor
299 300 301
      for ( size_t i = 0; i < valueCount; i++ )
        data[i] = (float)(data[i] * scaleFactor);
      break;
302 303 304 305 306 307 308 309 310 311 312
    case 1: // haveMissVal
      if (missValIsNaN)
        {
          for ( size_t i = 0; i < valueCount; i++ )
            missValCount += (size_t)DBL_IS_NAN(data[i]);
        }
      else
        {
          for ( size_t i = 0; i < valueCount; i++ )
            missValCount += (size_t)IS_EQUAL(data[i], missVal);
        }
313 314 315 316 317 318 319 320 321 322 323 324 325
      break;
    }

  return missValCount;
}

static
size_t min_size(size_t a, size_t b)
{
  return a < b ? a : b;
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
326
void transpose2dArrayDP(int gridId, double *data)
327
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328 329 330
  size_t inWidth = gridInqYsize(gridId);
  size_t inHeight = gridInqXsize(gridId);
  
331 332
  const size_t cacheBlockSize = 256;   // Purely an optimization parameter. Current value of 32 means we are handling 8kB blocks,
                                       // which should be a decent compromise on many architectures.
333 334
  double **out = (double**) malloc(inWidth*sizeof(double*));
  double **temp = (double**) malloc(inHeight*sizeof(double*));
335
  temp[0] = (double *) malloc(inHeight*inWidth*sizeof(double));
336
  memcpy(temp[0], data, inHeight*inWidth*sizeof(double));
337 338
  for (size_t i = 0; i < inWidth; i++) out[i] = data + (inHeight*i);
  for (size_t i = 1; i < inHeight; i++) temp[i] = temp[0] + (inWidth*i);
339 340 341 342 343 344 345 346

  /*
  for ( size_t y = 0; y < inHeight; ++y )
    for ( size_t x = 0; x < inWidth; ++x )
      out[x][y] = temp[y][x];
  */

  for ( size_t yBlock = 0; yBlock < inHeight; yBlock += cacheBlockSize )
347 348 349 350 351 352
    for ( size_t xBlock = 0; xBlock < inWidth; xBlock += cacheBlockSize )
      for ( size_t y = yBlock, yEnd = min_size(yBlock + cacheBlockSize, inHeight); y < yEnd; y++ )
        for ( size_t x = xBlock, xEnd = min_size(xBlock + cacheBlockSize, inWidth); x < xEnd; x++ )
          {
            out[x][y] = temp[y][x];
          }
353

354
  free(out);
355 356
  free(temp[0]);
  free(temp);
357 358 359
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360
void transpose2dArraySP(size_t gridId, float *data)
361
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
362 363 364
  size_t inWidth = gridInqYsize(gridId);
  size_t inHeight = gridInqXsize(gridId);

365 366
  const size_t cacheBlockSize = 256;   // Purely an optimization parameter. Current value of 32 means we are handling 8kB blocks,
                                       // which should be a decent compromise on many architectures.
367 368
  float **out = (float**) malloc(inWidth*sizeof(float*));
  float **temp = (float**) malloc(inHeight*sizeof(float*));
369
  temp[0] = (float *) malloc(inHeight*inWidth*sizeof(float));
370
  memcpy(temp[0], data, inHeight*inWidth*sizeof(float));
371 372
  for (size_t i = 0; i < inWidth; i++) out[i] = data + (inHeight*i);
  for (size_t i = 1; i < inHeight; i++) temp[i] = temp[0] + (inWidth*i);
373 374 375 376 377 378 379 380

  /*
  for ( size_t y = 0; y < inHeight; ++y )
    for ( size_t x = 0; x < inWidth; ++x )
      out[x][y] = temp[y][x];
  */

  for ( size_t yBlock = 0; yBlock < inHeight; yBlock += cacheBlockSize )
381 382 383 384 385 386
    for ( size_t xBlock = 0; xBlock < inWidth; xBlock += cacheBlockSize )
      for ( size_t y = yBlock, yEnd = min_size(yBlock + cacheBlockSize, inHeight); y < yEnd; y++ )
        for ( size_t x = xBlock, xEnd = min_size(xBlock + cacheBlockSize, inWidth); x < xEnd; x++ )
          {
            out[x][y] = temp[y][x];
          }
387

388
  free(out);
389 390
  free(temp[0]);
  free(temp);
391 392 393 394 395
}

static
void cdfInqDimIds(stream_t *streamptr, int varId, int (*outDimIds)[3])
{
396 397
  const int gridId = vlistInqVarGrid(streamptr->vlistID, varId);
  const int gridindex = vlistGridIndex(streamptr->vlistID, gridId);
398

399
  (*outDimIds)[0] = (*outDimIds)[1] = (*outDimIds)[2] = CDI_UNDEFID;
400 401 402 403 404 405 406
  switch ( gridInqType(gridId) )
    {
      case GRID_TRAJECTORY:
        cdfReadGridTraj(streamptr, gridId);
        break;

      case GRID_UNSTRUCTURED:
407
        (*outDimIds)[0] = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
408 409
        break;

410 411 412 413
      case GRID_GAUSSIAN_REDUCED:
        (*outDimIds)[0] = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
        break;

414
      default:
415
        (*outDimIds)[0] = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_X];
416
        (*outDimIds)[1] = streamptr->ncgrid[gridindex].ncIDs[CDF_DIMID_Y];
417 418 419
        break;
    }

420 421
  const int zaxisID = vlistInqVarZaxis(streamptr->vlistID, varId);
  const int zaxisindex = vlistZaxisIndex(streamptr->vlistID, zaxisID);
422 423 424 425 426 427
  (*outDimIds)[2] = streamptr->zaxisID[zaxisindex];
}

static
int cdfGetSkipDim(int fileId, int ncvarid, int (*dimIds)[3])
{
428 429
  if((*dimIds)[0] != CDI_UNDEFID) return 0;
  if((*dimIds)[1] != CDI_UNDEFID) return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
430

431 432 433 434 435 436
  int nvdims;
  cdf_inq_varndims(fileId, ncvarid, &nvdims);
  if(nvdims != 3) return 0;

  int varDimIds[3];
  cdf_inq_vardimid(fileId, ncvarid, varDimIds);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437

438 439 440 441 442 443 444 445 446 447 448
  size_t size = 0;
  if ( (*dimIds)[2] == varDimIds[2] )
    {
      cdf_inq_dimlen(fileId, varDimIds[1], &size);
      if ( size == 1 ) return 1;
    }
  else if ( (*dimIds)[2] == varDimIds[1] )
    {
      cdf_inq_dimlen(fileId, varDimIds[2], &size);
      if ( size == 1 ) return 2;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449

450 451 452 453
  return 0;
}

static
454
void cdfGetSliceSlapDescription(stream_t *streamptr, int varID, int levelID, bool *outSwapXY, size_t (*start)[4], size_t (*count)[4])
455
{
456
  const int tsID = streamptr->curTsID;
457 458
  if ( CDI_Debug ) Message("tsID = %d", tsID);

459
  const int fileId = streamptr->fileID;
460 461
  const int vlistID = streamptr->vlistID;
  const int ncvarid = streamptr->vars[varID].ncvarid;
462

463 464
  const int gridId = vlistInqVarGrid(vlistID, varID);
  const int timetype = vlistInqVarTimetype(vlistID, varID);
465
  const size_t gridsize = gridInqSize(gridId);
466 467 468

  streamptr->numvals += gridsize;

469
  int dimIds[3]; // this array joins the old variables xid, yid, and zid
470
  cdfInqDimIds(streamptr, varID, &dimIds);
471

472
  const int skipdim = cdfGetSkipDim(fileId, ncvarid, &dimIds);
473 474

  int dimorder[3];
475
  vlistInqVarDimorder(vlistID, varID, &dimorder);
476

477
  *outSwapXY = (dimorder[2] == 2 || dimorder[0] == 1) && dimIds[0] != CDI_UNDEFID && dimIds[1] != CDI_UNDEFID ;
478 479 480 481 482 483 484 485 486

  int ndims = 0;

#define addDimension(startIndex, extent) do {   \
      (*start)[ndims] = startIndex; \
      (*count)[ndims] = extent; \
      ndims++; \
  } while(0)

487
  if ( timetype != TIME_CONSTANT ) addDimension((size_t)tsID, 1);
488 489 490 491 492
  if ( skipdim == 1 ) addDimension(0, 1);

  for ( int id = 0; id < 3; ++id )
    {
      size_t size;
493
      const int curDimId = dimIds[dimorder[id]-1];
494
      if ( curDimId == CDI_UNDEFID ) continue;
495 496 497 498 499 500 501 502
      switch ( dimorder[id] )
        {
          case 1:
          case 2:
            cdf_inq_dimlen(fileId, curDimId, &size);
            addDimension(0, size);
            break;
          case 3:
503
            addDimension((size_t)levelID, 1);
504 505 506 507 508 509 510 511 512 513 514 515 516 517 518
            break;
          default:
            Error("Internal errror: Malformed dimension order encountered. Please report this bug.\n");
        }
    }

  if ( skipdim == 2 ) addDimension(0, 1);

  assert(ndims <= (int)(sizeof(*start)/sizeof(**start)));
  assert(ndims <= (int)(sizeof(*count)/sizeof(**count)));

#undef addDimension

  if ( CDI_Debug )
    for (int idim = 0; idim < ndims; idim++)
519
      Message("dim = %d  start = %d  count = %d", idim, (*start)[idim], (*count)[idim]);
520 521 522 523 524

  int nvdims;
  cdf_inq_varndims(fileId, ncvarid, &nvdims);

  if ( nvdims != ndims )
525 526 527 528 529
    {
      char name[CDI_MAX_NAME];
      vlistInqVarName(vlistID, varID, name);
      Error("Internal error, variable %s has an unsupported array structure!", name);
    }
530 531
}

532
static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
size_t getSizeVar3D(int vlistID, int varID)
534
{
535 536
  const int gridID  = vlistInqVarGrid(vlistID, varID);
  const int zaxisID = vlistInqVarZaxis(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587
  return gridInqSize(gridID) * (size_t)zaxisInqSize(zaxisID);
}


static
void cdfReadDataSliceSP2DP(int fileID, int ncvarid, size_t length, size_t start[4], size_t count[4], double *data)
{
  float *data_fp = (float *) Malloc(length*sizeof(*data_fp));
  cdf_get_vara_float(fileID, ncvarid, start, count, data_fp);
  for ( size_t i = 0; i < length; i++ ) data[i] = (double) data_fp[i];
  Free(data_fp);
}

static
void cdfReadDataSliceDP2SP(int fileID, int ncvarid, size_t length, size_t start[4], size_t count[4], float *data)
{
  double *data_dp = (double *) Malloc(length*sizeof(*data_dp));
  cdf_get_vara_double(fileID, ncvarid, start, count, data_dp);
  for ( size_t i = 0; i < length; i++ ) data[i] = (float) data_dp[i];
  Free(data_dp);
}

static
void cdfCheckDataDP_UINT8(int fileID, int ncvarid, int vlistID, int varID, size_t length, double *data)
{
  if ( vlistInqVarDatatype(vlistID, varID) == CDI_DATATYPE_UINT8 )
    {
      nc_type xtype;
      cdf_inq_vartype(fileID, ncvarid, &xtype);
      if ( xtype == NC_BYTE )
        {
          for ( size_t i = 0; i < length; i++ )
            if ( data[i] < 0 ) data[i] += 256;
        }
    }
}

static
void cdfCheckDataSP_UINT8(int fileID, int ncvarid, int vlistID, int varID, size_t length, float *data)
{
  if ( vlistInqVarDatatype(vlistID, varID) == CDI_DATATYPE_UINT8 )
    {
      nc_type xtype;
      cdf_inq_vartype(fileID, ncvarid, &xtype);
      if ( xtype == NC_BYTE )
        {
          for ( size_t i = 0; i < length; i++ )
            if ( data[i] < 0 ) data[i] += 256;
        }
    }
}
588

Uwe Schulzweida's avatar
Uwe Schulzweida committed
589 590 591
static
void cdfReadDataDP(stream_t *streamptr, int varID, size_t length, size_t start[4], size_t count[4], double *data)
{
592 593 594
  const int vlistID = streamptr->vlistID;
  const int fileID = streamptr->fileID;
  const int ncvarid = streamptr->vars[varID].ncvarid;
595
  const int datatype = vlistInqVarDatatype(vlistID, varID);
596

597
  if (datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
598
    {
599 600 601 602 603 604 605 606 607
      cdf_get_vara(fileID, ncvarid, start, count, data);
      if (datatype == CDI_DATATYPE_CPX32)
        {
          for (long i = (long)length-1; i >= 0; i--)
            {
              data[2*i] = (double)(((float*)data)[2*i]);
              data[2*i+1] = (double)(((float*)data)[2*i+1]);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
608 609 610
    }
  else
    {
611 612 613 614 615 616 617
      if ( datatype == CDI_DATATYPE_FLT32 )
        {
          cdfReadDataSliceSP2DP(fileID, ncvarid, length, start, count, data);
        }
      else
        {
          cdf_get_vara_double(fileID, ncvarid, start, count, data);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
      
619 620
          cdfCheckDataDP_UINT8(fileID, ncvarid, vlistID, varID, length, data);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
621 622 623 624 625 626
    }
}

static
void cdfReadDataSP(stream_t *streamptr, int varID, size_t length, size_t start[4], size_t count[4], float *data)
{
627 628 629
  const int vlistID = streamptr->vlistID;
  const int fileID = streamptr->fileID;
  const int ncvarid = streamptr->vars[varID].ncvarid;
630
  const int datatype = vlistInqVarDatatype(vlistID, varID);
631

632
  if (datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
    {
634 635 636 637 638 639 640 641 642 643 644 645 646 647 648
      if (datatype == CDI_DATATYPE_CPX64)
        {
          double *cdata = (double *) Malloc(2*length*sizeof(double));
          cdf_get_vara(fileID, ncvarid, start, count, cdata);
          for (size_t i = 0; i < length; i++)
            {
              data[2*i] = (float)(cdata[2*i]);
              data[2*i+1] = (float)(cdata[2*i+1]);
            }
          Free(cdata);
        }
      else
        {
          cdf_get_vara(fileID, ncvarid, start, count, data);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
649 650 651
    }
  else
    {
652 653 654 655 656 657 658
      if ( datatype == CDI_DATATYPE_FLT64 )
        {
          cdfReadDataSliceDP2SP(fileID, ncvarid, length, start, count, data);
        }
      else
        {
          cdf_get_vara_float(fileID, ncvarid, start, count, data);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
659
      
660 661
          cdfCheckDataSP_UINT8(fileID, ncvarid, vlistID, varID, length, data);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
662 663 664 665 666 667 668 669
    }
}

static
void cdfReadVarDP(stream_t *streamptr, int varID, double *data, size_t *nmiss)
{
  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamptr->self, varID);

670
  const int vlistID = streamptr->vlistID;
671

Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
  size_t start[4], count[4];
673 674
  cdfGetSlapDescription(streamptr, varID, &start, &count);

675
  const size_t length = getSizeVar3D(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
  cdfReadDataDP(streamptr, varID, length, start, count, data);
677

Uwe Schulzweida's avatar
Uwe Schulzweida committed
678
  *nmiss = cdfDoInputDataTransformationDP(vlistID, varID, length, data);
679 680
}

681
static
682
void cdfReadVarSP(stream_t *streamptr, int varID, float *data, size_t *nmiss)
683 684 685
{
  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamptr->self, varID);

686
  const int vlistID = streamptr->vlistID;
687

Uwe Schulzweida's avatar
Uwe Schulzweida committed
688
  size_t start[4], count[4];
689 690
  cdfGetSlapDescription(streamptr, varID, &start, &count);

691
  const size_t length = getSizeVar3D(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
692
  cdfReadDataSP(streamptr, varID, length, start, count, data);
693

Uwe Schulzweida's avatar
Uwe Schulzweida committed
694
  *nmiss = cdfDoInputDataTransformationSP(vlistID, varID, length, data);
695 696 697
}


698
void cdf_read_var(stream_t *streamptr, int varID, int memtype, void *data, size_t *nmiss)
699 700 701 702 703 704 705 706
{
  if ( memtype == MEMTYPE_DOUBLE )
    cdfReadVarDP(streamptr, varID, (double*) data, nmiss);
  else
    cdfReadVarSP(streamptr, varID, (float*) data, nmiss);
}

static
707
void cdfReadVarSliceDP(stream_t *streamptr, int varID, int levelID, double *data, size_t *nmiss)
708 709 710 711 712
{
  if ( CDI_Debug )
    Message("streamID = %d  varID = %d  levelID = %d", streamptr->self, varID, levelID);

  bool swapxy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
713
  size_t start[4], count[4];
714 715
  cdfGetSliceSlapDescription(streamptr, varID, levelID, &swapxy, &start, &count);

716 717 718
  const int vlistID = streamptr->vlistID;
  const int gridId = vlistInqVarGrid(vlistID, varID);
  const size_t length = gridInqSize(gridId);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
719
  cdfReadDataDP(streamptr, varID, length, start, count, data);
720

Uwe Schulzweida's avatar
Uwe Schulzweida committed
721
  if ( swapxy ) transpose2dArrayDP(gridId, data);
722

Uwe Schulzweida's avatar
Uwe Schulzweida committed
723
  *nmiss = cdfDoInputDataTransformationDP(vlistID, varID, length, data);
724 725
}

726
static
727
void cdfReadVarSliceSP(stream_t *streamptr, int varID, int levelID, float *data, size_t *nmiss)
728 729 730 731 732
{
  if ( CDI_Debug )
    Message("streamID = %d  varID = %d  levelID = %d", streamptr->self, varID, levelID);

  bool swapxy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
  size_t start[4], count[4];
734 735
  cdfGetSliceSlapDescription(streamptr, varID, levelID, &swapxy, &start, &count);

736 737 738
  const int vlistID = streamptr->vlistID;
  const int gridId = vlistInqVarGrid(vlistID, varID);
  const size_t length = gridInqSize(gridId);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
739
  cdfReadDataSP(streamptr, varID, length, start, count, data);
740

Uwe Schulzweida's avatar
Uwe Schulzweida committed
741
  if ( swapxy ) transpose2dArraySP(gridId, data);
742

Uwe Schulzweida's avatar
Uwe Schulzweida committed
743
  *nmiss = cdfDoInputDataTransformationSP(vlistID, varID, length, data);
744 745 746
}


747
void cdf_read_var_slice(stream_t *streamptr, int varID, int levelID, int memtype, void *data, size_t *nmiss)
748 749 750 751 752 753 754 755
{
  if ( memtype == MEMTYPE_DOUBLE )
    cdfReadVarSliceDP(streamptr, varID, levelID, (double*) data, nmiss);
  else
    cdfReadVarSliceSP(streamptr, varID, levelID, (float*) data, nmiss);
}


756
void cdf_read_record(stream_t *streamptr, int memtype, void *data, size_t *nmiss)
757 758 759
{
  if ( CDI_Debug ) Message("streamID = %d", streamptr->self);

760 761 762 763 764
  const int tsID    = streamptr->curTsID;
  const int vrecID  = streamptr->tsteps[tsID].curRecID;
  const int recID   = streamptr->tsteps[tsID].recIDs[vrecID];
  const int varID   = streamptr->tsteps[tsID].records[recID].varID;
  const int levelID = streamptr->tsteps[tsID].records[recID].levelID;
765

766
  if ( memtype == MEMTYPE_DOUBLE )
767
    cdfReadVarSliceDP(streamptr, varID, levelID, (double*) data, nmiss);
768
  else
769
    cdfReadVarSliceSP(streamptr, varID, levelID, (float*) data, nmiss);
770 771
}

772 773 774 775 776 777 778 779 780 781 782
//----------------------------------------------------------------------------
// Parallel Version
//----------------------------------------------------------------------------

void cdfReadVarSliceDPPart(stream_t *streamptr, int varID, int levelID, int varType, int startpoint, size_t length, double *data, size_t *nmiss)
{
  (void)(varType);

  if ( CDI_Debug )
    Message("streamID = %d  varID = %d  levelID = %d", streamptr->self, varID, levelID);

783
  const int vlistID = streamptr->vlistID;
784 785

  bool swapxy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
786
  size_t start[4], count[4];
787 788
  cdfGetSliceSlapDescription(streamptr, varID, levelID, &swapxy, &start, &count);

789 790
  const int gridId = vlistInqVarGrid(vlistID, varID);
  const size_t gridsize = gridInqSize(gridId);
791 792

  unsigned int position = 0; 
793 794
  for (int i = 0; i < 4; i++)
    if (count[i] == gridsize) position = i;
795 796 797 798

  start[position] = start[position]+startpoint;
  count[position] = length;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
799
  cdfReadDataDP(streamptr, varID, length, start, count, data);
800

Uwe Schulzweida's avatar
Uwe Schulzweida committed
801
  if ( swapxy ) transpose2dArrayDP(gridId, data);
802

Uwe Schulzweida's avatar
Uwe Schulzweida committed
803
  *nmiss = cdfDoInputDataTransformationDP(vlistID, varID, length, data);
804 805 806 807 808 809 810 811 812
}

void cdfReadVarSliceSPPart(stream_t *streamptr, int varID, int levelID, int varType, int startpoint, size_t length, float *data, size_t *nmiss)
{
  (void)(varType);

  if ( CDI_Debug )
    Message("streamID = %d  varID = %d  levelID = %d", streamptr->self, varID, levelID);

813
  const int vlistID = streamptr->vlistID;
814 815

  bool swapxy;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
816
  size_t start[4], count[4];
817 818
  cdfGetSliceSlapDescription(streamptr, varID, levelID, &swapxy, &start, &count);

819
  const int gridId = vlistInqVarGrid(vlistID, varID);
820
  size_t gridsize = gridInqSize(gridId);
821 822

  unsigned int position = 0; 
823 824
  for (int i = 0; i < 4; i++)
    if (count[i] == gridsize) position = i;
825 826 827 828

  start[position] = start[position]+startpoint;
  count[position] = length;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
829
  cdfReadDataSP(streamptr, varID, length, start, count, data);
830

Uwe Schulzweida's avatar
Uwe Schulzweida committed
831
  if ( swapxy ) transpose2dArraySP(gridId, data);
832

Uwe Schulzweida's avatar
Uwe Schulzweida committed
833
  *nmiss = cdfDoInputDataTransformationSP(vlistID, varID, length, data);
834 835 836 837 838 839 840 841 842 843 844 845 846
}

static
int cdiStreamReadVarSlicePart(int streamID, int varID, int levelID, int varType, int start, size_t size, int memtype, void *data, size_t *nmiss)
{
  int status = 0;

  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamID, varID);

  check_parg(data);
  check_parg(nmiss);

  stream_t *streamptr = stream_to_pointer(streamID);
847
  const int filetype = streamptr->filetype;
848 849 850 851 852 853

  *nmiss = 0;

  // currently we only care for netcdf data
  switch (filetype)
    {
854 855 856 857 858 859 860 861 862
#ifdef HAVE_LIBGRIB
    case CDI_FILETYPE_GRB:
    case CDI_FILETYPE_GRB2:
      {
        grb_read_var_slice(streamptr, varID, levelID, memtype, data, nmiss);
	break;
      }
#endif
#ifdef HAVE_LIBNETCDF
863 864 865 866 867
    case CDI_FILETYPE_NC:
    case CDI_FILETYPE_NC2:
    case CDI_FILETYPE_NC4:
    case CDI_FILETYPE_NC4C:
    case CDI_FILETYPE_NC5:
868 869
      {
        if ( memtype == MEMTYPE_FLOAT )
870
          cdfReadVarSliceSPPart(streamptr, varID, levelID, varType, start, size, (float *)data, nmiss);
871
        else
872
          cdfReadVarSliceDPPart(streamptr, varID, levelID, varType, start, size, (double *)data, nmiss);
873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892
        break;
      }
#endif
    default:
      {
        Error("%s support not compiled in!", strfiletype(filetype));
        status = 2;
        break;
      }
    }

  return status;
}


void cdfReadVarDPPart(stream_t *streamptr, int varID, int varType, int startpoint, size_t length, double *data, size_t *nmiss)
{
  (void)(varType);
  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamptr->self, varID);

893 894
  const int vlistID = streamptr->vlistID;
  const int ncvarid = streamptr->vars[varID].ncvarid;
895

Uwe Schulzweida's avatar
Uwe Schulzweida committed
896
  size_t start[4], count[4];
897 898
  cdfGetSlapDescription(streamptr, varID, &start, &count);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
899 900 901
  const int ltime = TIME_CONSTANT != vlistInqVarTimetype(vlistID, varID);
  start[1+ltime] = start[1+ltime]+startpoint;
  count[1+ltime] = length;
902

Uwe Schulzweida's avatar
Uwe Schulzweida committed
903
  cdf_get_vara_double(streamptr->fileID, ncvarid, start, count, data);
904

Uwe Schulzweida's avatar
Uwe Schulzweida committed
905
  *nmiss = cdfDoInputDataTransformationDP(vlistID, varID, length, data);
906 907 908 909 910 911 912
}

void cdfReadVarSPPart(stream_t *streamptr, int varID, int varType, int startpoint, size_t length, float *data, size_t *nmiss)
{
  (void)(varType);
  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamptr->self, varID);

913 914
  const int vlistID = streamptr->vlistID;
  const int ncvarid = streamptr->vars[varID].ncvarid;
915

Uwe Schulzweida's avatar
Uwe Schulzweida committed
916
  size_t start[4], count[4];
917 918
  cdfGetSlapDescription(streamptr, varID, &start, &count);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
919 920 921
  const int ltime = TIME_CONSTANT != vlistInqVarTimetype(vlistID, varID);
  start[1+ltime] = start[1+ltime]+startpoint;
  count[1+ltime] = length;
922

Uwe Schulzweida's avatar
Uwe Schulzweida committed
923
  cdf_get_vara_float(streamptr->fileID, ncvarid, start, count, data);
924

Uwe Schulzweida's avatar
Uwe Schulzweida committed
925
  *nmiss = cdfDoInputDataTransformationSP(vlistID, varID, length, data);
926 927 928 929 930 931 932 933 934 935 936 937
}

static
void cdiStreamReadVarPart(int streamID, int varID, int varType, int start, size_t size, int memtype, void *data, size_t *nmiss)
{
  (void)(varType);
  if ( CDI_Debug ) Message("streamID = %d  varID = %d", streamID, varID);

  check_parg(data);
  check_parg(nmiss);

  stream_t *streamptr = stream_to_pointer(streamID);
938
  const int filetype = streamptr->filetype;
939 940 941 942 943 944

  *nmiss = 0;

  // currently we only care for netcdf data
  switch (filetype)
    {
945 946 947 948 949 950 951 952 953
#ifdef HAVE_LIBGRIB
    case CDI_FILETYPE_GRB:
    case CDI_FILETYPE_GRB2:
      {
        grb_read_var(streamptr, varID, memtype, data, nmiss);
	break;
      }
#endif
#ifdef HAVE_LIBNETCDF
954 955 956 957 958
    case CDI_FILETYPE_NC:
    case CDI_FILETYPE_NC2:
    case CDI_FILETYPE_NC4:
    case CDI_FILETYPE_NC4C:
    case CDI_FILETYPE_NC5:
959 960
      {
        if ( memtype == MEMTYPE_FLOAT )
961
          cdfReadVarSPPart(streamptr, varID, varType, start, size, (float *)data, nmiss);
962
        else
963
          cdfReadVarDPPart(streamptr, varID, varType, start, size, (double *)data, nmiss);
964 965 966 967 968 969 970 971 972 973 974 975 97