Gradsdes.cc

/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2018 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
  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.
*/

/*
   This module contains the following operators:

      Gradsdes   gradsdes       GrADS data descriptor file
*/

#include <cdi.h>

#include "cdo_int.h"
#include "grid.h"
#include "pstream_int.h"
#include "util_fileextensions.h"

/*
  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

*/
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  */
};
struct gaindxb
{
  int bignum;    /* Number of off_t values */
  off_t *bigpnt; /* Pointer to off_t values */
};

/* Byte swap requested number of 4 byte elements */

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;
    }
}

/*
 * convert an IBM float to single precision number v1.0
 *
 *                      Wesley Ebisuzaki
 */
static float
ibm2flt(unsigned char *ibm)
{

  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;
  while (abspower)
    {
      if (abspower & 1)
        {
          value *= exp;
        }
      exp = exp * exp;
      abspower >>= 1;
    }

  if (power < 0) value = 1.0 / value;
  value = value * mant / 16777216.0;
  if (positive == 0) value = -value;
  return (float) value;
}

/*
 * convert a float to an IBM single precision number v1.0
 *
 *                      Wesley Ebisuzaki
 *
 * doesn't handle subnormal numbers
 */
static int
flt2ibm(float x, unsigned char *ibm)
{

  int sign, exp, i;
  double mant;

  if (!(fabs((double) x) > 0))
    {
      ibm[0] = ibm[1] = ibm[2] = ibm[3] = 0;
      return 0;
    }

  /* sign bit */
  if (x < 0.0)
    {
      sign = 128;
      x = -x;
    }
  else
    sign = 0;

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

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

  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;
    }

  /* 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;

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

  return 0;
}

#define GET_UINT4(a, b, c, d) ((int) ((a << 24) + (b << 16) + (c << 8) + (d)))
#define Put1Byte(buf, cnt, ival) (buf[cnt++] = (ival))
#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)))

#define PutInt(buf, cnt, ival) (ival < 0 ? Put4Byte(buf, cnt, 0x7fffffff - ival + 1) : Put4Byte(buf, cnt, ival))

static void
dumpmap()
{
  unsigned char urec[4];
  unsigned char vermap;
  unsigned char mrec[512];
  int swpflg = 0;
  int i;
  int nrecords = 0;
  struct gaindx indx;
  struct gaindxb indxb;
  size_t nbytes;
  FILE *mapfp;

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

  mapfp = fopen(cdoGetStreamName(0), "r");
  if (mapfp == NULL) cdoAbort("Open failed on %s", cdoGetStreamName(0));

  /* check the version number */

  fseek(mapfp, 1, 0);
  nbytes = fread(&vermap, sizeof(unsigned char), 1, mapfp);

  if (vermap == 0) vermap = 1;

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

  if (vermap == 2)
    {
      fseek(mapfp, 2, 0);

      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
      indx.hinum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);

      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
      indx.hfnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);

      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
      indx.intnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);

      nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
      indx.fltnum = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);

      nbytes = fread(mrec, sizeof(unsigned char), 7, mapfp);

      if (indx.hinum > 0)
        {
          indx.hipnt = (int *) Malloc(sizeof(int) * indx.hinum);
          for (i = 0; i < indx.hinum; i++)
            {
              nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
              indx.hipnt[i] = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
            }
        }
      if (indx.hfnum > 0)
        {
          indx.hfpnt = (float *) Malloc(sizeof(float) * indx.hfnum);
          nbytes = fread(indx.hfpnt, sizeof(float), indx.hfnum, mapfp);
        }
      if (indx.intnum > 0)
        {
          indx.intpnt = (int *) Malloc(sizeof(int) * indx.intnum);
          for (i = 0; i < indx.intnum; i++)
            {
              nbytes = fread(mrec, sizeof(unsigned char), 4, mapfp);
              indx.intpnt[i] = GET_UINT4(mrec[0], mrec[1], mrec[2], mrec[3]);
              if (indx.intpnt[i] < 0) indx.intpnt[i] = 0x7fffffff - indx.intpnt[i] + 1;
            }
        }
      if (indx.fltnum > 0)
        {
          indx.fltpnt = (float *) Malloc(sizeof(float) * indx.fltnum);
          for (i = 0; i < indx.fltnum; i++)
            {
              nbytes = fread(urec, sizeof(unsigned char), 4, mapfp);
              indx.fltpnt[i] = ibm2flt(urec);
            }
        }
    }
  else
    {
      fseek(mapfp, 0, 0);
      nbytes = fread(&indx, sizeof(struct gaindx), 1, mapfp);

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

      if (indx.hinum > 0)
        {
          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);
        }
      if (indx.hfnum > 0)
        {
          indx.hfpnt = (float *) Malloc(sizeof(float) * indx.hfnum);
          nbytes = fread(indx.hfpnt, sizeof(float), indx.hfnum, mapfp);
          if (swpflg) gabswp(indx.hfpnt, indx.hfnum);
        }

      if (indx.intnum > 0)
        {
          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);
        }
      if (indx.fltnum > 0)
        {
          indx.fltpnt = (float *) Malloc(sizeof(float) * indx.fltnum);
          nbytes = fread(indx.fltpnt, sizeof(float), indx.fltnum, mapfp);
          if (swpflg) gabswp(indx.fltpnt, indx.fltnum);
        }

      if (indx.hipnt[0] == 4)
        {
          indxb.bignum = indx.hipnt[4];
          if (indxb.bignum > 0)
            {
              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);
            }
        }
    }

  UNUSED(nbytes);

  fclose(mapfp);

  printf("hinum: %d\n", indx.hinum);
  for (i = 0; i < indx.hinum; i++) printf("%3d %5d\n", i + 1, indx.hipnt[i]);

  printf("\n");
  printf("hfnum: %d\n", indx.hfnum);
  for (i = 0; i < indx.hfnum; i++) printf("%3d %g\n", i + 1, indx.hfpnt[i]);

  printf("\n");

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

  if (indx.intnum == indx.fltnum)
    {
      printf("num: %d\n", indx.intnum);
      for (i = 0; i < indx.intnum / 3; i++)
        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]);
    }
  else if (indx.intnum == nrecords && indx.fltnum == nrecords * 3 && indxb.bignum == nrecords * 2)
    {
      printf("nrecords: %d\n", nrecords);
      for (i = 0; i < nrecords; i++)
        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]);
    }
  else
    {
      printf("intnum: %d\n", indx.intnum);
      for (i = 0; i < indx.intnum; i++) printf("%3d %d\n", i + 1, indx.intpnt[i]);

      printf("\n");
      printf("fltnum: %d\n", indx.fltnum);
      for (i = 0; i < indx.fltnum; i++) printf("%3d %g\n", i + 1, indx.fltpnt[i]);

      printf("\n");
      printf("bignum: %d\n", indxb.bignum);
      for (i = 0; i < indxb.bignum; i++) printf("%3d %zd\n", i + 1, (size_t) indxb.bigpnt[i]);
    }
}

static void
ctl_xydef(FILE *gdp, int gridID, bool *yrev)
{
  int i, j;
  double xfirst, yfirst, xinc, yinc;

  *yrev = false;

  int xsize = gridInqXsize(gridID);
  int ysize = gridInqYsize(gridID);

  int gridtype = gridInqType(gridID);
  int projtype = gridInqProjType(gridID);

  /* XDEF */

  if (gridtype == GRID_PROJECTION && projtype == CDI_PROJ_LCC)
    {
      double xmin = 1.e10, xmax = -1.e10, ymin = 1.e10, ymax = -1.e10;
      double xrange, yrange;
      int nx, ny, ni;
      double inc[] = { 1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001 };

      double xinc = gridInqXinc(gridID);
      double yinc = gridInqYinc(gridID);

      double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;

      gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
      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);

      gridID = gridToCurvilinear(gridID, 0);
      std::vector<double> xvals(xsize * ysize);
      std::vector<double> yvals(xsize * ysize);
      gridInqXvals(gridID, xvals.data());
      gridInqYvals(gridID, yvals.data());
      for (i = 0; i < xsize * ysize; ++i)
        {
          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];
        }

      xfirst = ((int) (xmin - 0.0));
      yfirst = ((int) (ymin - 0.0));
      xrange = ((int) (xmax + 1.5)) - xfirst;
      yrange = ((int) (ymax + 1.5)) - yfirst;

      ni = sizeof(inc) / sizeof(inc[0]);
      for (i = 0; i < ni; i++)
        {
          xinc = yinc = inc[i];
          nx = 1 + (int) (xrange / xinc);
          ny = 1 + (int) (yrange / yinc);

          if (nx > 1.5 * xsize && ny > 1.5 * ysize) break;
        }

      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);
      xinc = gridInqXinc(gridID);
      if (IS_EQUAL(xinc, 0) && gridInqXvals(gridID, NULL))
        {
          std::vector<double> xvals(xsize);
          gridInqXvals(gridID, xvals.data());
          fprintf(gdp, "XDEF %d LEVELS ", xsize);
          j = 0;
          for (i = 0; i < xsize; i++)
            {
              fprintf(gdp, "%7.3f ", xvals[i]);
              j++;
              if (j == 6)
                {
                  fprintf(gdp, "\n");
                  j = 0;
                  if (i != xsize - 1) fprintf(gdp, "               ");
                }
            }
          if (j) fprintf(gdp, "\n");
        }
      else
        {
          if (IS_EQUAL(xinc, 0)) xinc = 360.0 / xsize;
          fprintf(gdp, "XDEF %d LINEAR %f %f\n", xsize, xfirst, xinc);
        }
    }

  /* YDEF */

  if (!(gridtype == GRID_PROJECTION && projtype == CDI_PROJ_LCC))
    {
      yfirst = gridInqYval(gridID, 0);
      yinc = gridInqYinc(gridID);
      if (gridtype == GRID_GAUSSIAN) yinc = 0;

      if (IS_EQUAL(yinc, 0) && gridInqYvals(gridID, NULL))
        {
          std::vector<double> yvals(ysize);
          gridInqYvals(gridID, yvals.data());
          fprintf(gdp, "YDEF %d LEVELS ", ysize);
          j = 0;
          if (yvals[0] > yvals[ysize - 1])
            {
              *yrev = true;
              for (i = ysize - 1; i >= 0; i--)
                {
                  fprintf(gdp, "%7.3f ", yvals[i]);
                  j++;
                  if (j == 6)
                    {
                      fprintf(gdp, "\n");
                      j = 0;
                      if (i != 0) fprintf(gdp, "               ");
                    }
                }
            }
          else
            {
              for (i = 0; i < ysize; i++)
                {
                  fprintf(gdp, "%7.3f ", yvals[i]);
                  j++;
                  if (j == 6)
                    {
                      fprintf(gdp, "\n");
                      j = 0;
                      if (i != ysize - 1) fprintf(gdp, "               ");
                    }
                }
            }

          if (j) fprintf(gdp, "\n");
        }
      else
        {
          if (IS_EQUAL(yinc, 0)) yinc = 180.0 / ysize;
          if (yinc < 0)
            {
              *yrev = true;
              fprintf(gdp, "YDEF %d LINEAR %f %f\n", ysize, yfirst + yinc * (ysize - 1), -yinc);
            }
          else
            fprintf(gdp, "YDEF %d LINEAR %f %f\n", ysize, yfirst, yinc);
        }
    }
}

static void
ctl_zdef(FILE *gdp, int vlistID, bool *zrev)
{
  int i, j, index;
  int zaxisIDmax = -1;
  int zaxisID, nlev;
  double levinc = 0;

  *zrev = false;
  int nzaxis = vlistNzaxis(vlistID);

  int nlevmax = 0;
  for (index = 0; index < nzaxis; index++)
    {
      zaxisID = vlistZaxis(vlistID, index);
      nlev = zaxisInqSize(zaxisID);
      if (nlev > nlevmax)
        {
          nlevmax = nlev;
          zaxisIDmax = zaxisID;
        }
    }

  std::vector<double> levels(nlevmax);
  cdoZaxisInqLevels(zaxisIDmax, levels.data());
  bool lplev = (zaxisInqType(zaxisIDmax) == ZAXIS_PRESSURE);
  double level0 = levels[0];
  if (nlevmax > 1)
    {
      if (levels[0] < levels[1] && zaxisInqType(zaxisIDmax) != ZAXIS_HYBRID) *zrev = true;

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

      if (IS_EQUAL(levinc, 1)) *zrev = false;

      for (i = 1; i < nlevmax; i++)
        {
          if (IS_NOT_EQUAL(levinc, (levels[i] - levels[i - 1])))
            {
              levinc = 0;
              break;
            }
        }
    }

  if (IS_NOT_EQUAL(levinc, 0))
    fprintf(gdp, "ZDEF %d LINEAR %g %g\n", nlevmax, level0, levinc);
  else
    {
      fprintf(gdp, "ZDEF %d LEVELS ", nlevmax);
      j = 0;
      /* zrev not needed !!!
      if ( *zrev )
        {
          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, "               ");
                }
            }
        }
      else
      */
      {
        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");
    }
}

static void
ctl_options(FILE *gdp, bool yrev, bool zrev, bool sequential, bool bigendian, bool littleendian, bool flt64, bool cal365day)
{
  /* if ( filetype == CDI_FILETYPE_GRB ) zrev = false; */

  if (yrev || zrev || sequential || bigendian || littleendian || flt64)
    {
      fprintf(gdp, "OPTIONS");
      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");
      fprintf(gdp, "\n");
    }
}

static void
ctl_undef(FILE *gdp, int vlistID)
{
  double missval = vlistInqVarMissval(vlistID, 0);
  fprintf(gdp, "UNDEF  %g\n", missval);
}

static void
ctl_vars(FILE *gdp, int filetype, int vlistID, int nvarsout, int *vars)
{
  int ltype, code;
  int zaxisID, nlev;
  int i, j;
  int len;
  char varname[CDI_MAX_NAME], varlongname[CDI_MAX_NAME], varunits[CDI_MAX_NAME];

  int nvars = vlistNvars(vlistID);

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

  for (int varID = 0; varID < nvars; varID++)
    {
      if (vars[varID] == TRUE)
        {
          zaxisID = vlistInqVarZaxis(vlistID, varID);
          ltype = zaxisInqLtype(zaxisID);
          nlev = zaxisInqSize(zaxisID);
          vlistInqVarName(vlistID, varID, varname);

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

          if (i < len)
            for (j = i; j < len; j++) varname[j] = varname[j + 1];

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

          if (nlev == 1) nlev = 0;

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

          if (filetype == CDI_FILETYPE_GRB)
            {
              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);
            }
          else if (filetype == CDI_FILETYPE_NC)
            {
              int xyz = vlistInqVarXYZ(vlistID, varID);

              fprintf(gdp, "  ");
              if (vlistInqVarTimetype(vlistID, varID) != TIME_CONSTANT) fprintf(gdp, "t,");
              if (xyz == 321)
                {
                  if (nlev > 0) fprintf(gdp, "z,");
                  fprintf(gdp, "y,x");
                }
              else if (xyz == 312)
                {
                  if (nlev > 0) fprintf(gdp, "z,");
                  fprintf(gdp, "x,y");
                }
              else if (xyz == 231)
                {
                  fprintf(gdp, "y,");
                  if (nlev > 0) fprintf(gdp, "z,");
                  fprintf(gdp, "x");
                }
              else if (xyz == 132)
                {
                  fprintf(gdp, "x,");
                  if (nlev > 0) fprintf(gdp, "z,");
                  fprintf(gdp, "y");
                }
              else
                {
                  if (nlev > 0) fprintf(gdp, "z,");
                  fprintf(gdp, "y,x");
                }
            }
          else
            fprintf(gdp, "  99");

          if (varlongname[0] == 0)
            fprintf(gdp, "  %s", varname);
          else
            fprintf(gdp, "  %s", varlongname);

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

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

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

static void
write_map_grib1(const char *ctlfile, int map_version, int nrecords, int *intnum, float *fltnum, off_t *bignum)
{
  int i;
  struct gaindx indx;
  struct gaindxb indxb;
  FILE *mapfp;
  int hinum[5];

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

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

  indx.type = map_version;
  indx.hfnum = 0;
  if (map_version == 4)
    {
      indx.hinum = 5;
      indx.intnum = nrecords;
      indxb.bignum = 2 * nrecords;
    }
  else
    {
      indx.hinum = 4;
      indx.intnum = 3 * nrecords;
      indxb.bignum = 0;
    }
  indx.fltnum = 3 * nrecords;

  indx.hipnt = NULL;
  indx.hfpnt = NULL;
  indx.intpnt = NULL;
  indx.fltpnt = NULL;
  indxb.bigpnt = NULL;

  hinum[0] = map_version;
  hinum[1] = 1;
  hinum[2] = nrecords;
  hinum[3] = 255;
  hinum[4] = indxb.bignum;

  if (map_version == 2)
    {
      int nb, bcnt, rc, j;
      float fdum;
      unsigned char ibmfloat[4];

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

      nb = 2 + (indx.hinum * 4) + /* version in byte 2, then 4 ints with number of each data type */
           indx.hinum * sizeof(int) + indx.hfnum * sizeof(int) + indx.intnum * sizeof(int) + indx.fltnum * sizeof(float);

      /* add additional info */

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

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

      bcnt = 0;
      Put1Byte(map, bcnt, 0);
      Put1Byte(map, bcnt, map_version);

      Put4Byte(map, bcnt, indx.hinum);
      Put4Byte(map, bcnt, indx.hfnum);
      Put4Byte(map, bcnt, indx.intnum);
      Put4Byte(map, bcnt, indx.fltnum);

      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 */

      if (indx.hinum)
        for (i = 0; i < indx.hinum; i++) Put4Byte(map, bcnt, hinum[i]);

      if (indx.hfnum)
        {
          /* blank for now */
        }

      for (i = 0; i < nrecords; i++)
        {
          PutInt(map, bcnt, (int) bignum[i * 2]);
          PutInt(map, bcnt, (int) bignum[i * 2 + 1]);
          PutInt(map, bcnt, intnum[i]);
        }

      for (i = 0; i < indx.fltnum; i++)
        {
          fdum = fltnum[i];
          rc = flt2ibm(fdum, ibmfloat);
          if (rc < 0) cdoAbort("overflow in IBM float conversion");
          for (j = 0; j < 4; j++) map[bcnt++] = ibmfloat[j];
        }

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

      for (i = 0; i < 8; i++)
        {
          fdum = 0;
          rc = flt2ibm(fdum, ibmfloat);
          if (rc < 0) cdoAbort("overflow in IBM float conversion");
          for (j = 0; j < 4; j++) map[bcnt++] = ibmfloat[j];
        }

      fwrite(map, 1, bcnt, mapfp);

      Free(map);
    }
  else
    {
      fwrite(&indx, sizeof(struct gaindx), 1, mapfp);
      if (indx.hinum > 0) fwrite(hinum, sizeof(int), indx.hinum, mapfp);
      if (map_version == 1)
        {
          std::vector<int> intnumbuf(indx.intnum);
          for (i = 0; i < nrecords; i++)
            {
              intnumbuf[i * 3 + 0] = (int) bignum[i * 2];
              intnumbuf[i * 3 + 1] = (int) bignum[i * 2 + 1];
              intnumbuf[i * 3 + 2] = intnum[i];
            }
          if (indx.intnum > 0) fwrite(intnumbuf.data(), sizeof(int), indx.intnum, mapfp);
          if (indx.fltnum > 0) fwrite(fltnum, sizeof(float), indx.fltnum, mapfp);
        }
      else
        {
          if (indx.intnum > 0) fwrite(intnum, sizeof(int), indx.intnum, mapfp);
          if (indx.fltnum > 0) fwrite(fltnum, sizeof(float), indx.fltnum, mapfp);
          if (indxb.bignum > 0) fwrite(bignum, sizeof(off_t), indxb.bignum, mapfp);
        }
    }

  fclose(mapfp);
}

/*
static
void write_map_grib2(const char *ctlfile, int map_version, int nrecords, int
*intnum, float *fltnum, off_t *bignum)
{
}
*/

void *
Gradsdes(void *process)
{
  int gridID