/*
** File:  packgrib_.c
**
** Author:  Bob Dattore
**          NCAR/DSS
**          dattore@ucar.edu
**          (303) 497-1825
**
** Latest Revision:   2 Mar 99
**
** Purpose: to provide a single routine for packing GRIB grids that is callable
**          from Fortran
**
** Notes:   1) There are several routines defined in this file, but the only one
**             that is callable from Fortran is PACKGRIB.
**
**          2) The user is expected to have some understanding of the GRIB
**             format, as some of the input to PACKGRIB consists of codes, etc.
**             that are defined by the GRIB format.  You can find a description
**             of the GRIB format at http://www.scd.ucar.edu/dss/docs/gribdoc,
**             or you can get it by anonymous FTP at ncardata.ucar.edu, under
**             "datasets/ds083.0".
**
**          3) PACKGRIB has been tested on a Sun running Solaris and an SGI
**             running IRIX - please report any problems to dattore@ucar.edu.
**
**
** Fortran calling syntax:
**    INTEGER  GRIB_LBL(39)
**    INTEGER  PDS_EXT(IPDIM)
**    REAL     GRIDPOINTS(NUM_I,NUM_J)
**    INTEGER  OUTPUT_BUFFER(ODIM)
** 
**    CALL PACKGRIB(GRIB_LBL,PDS_EXT,IPDIM,GRIDPOINTS,IGIDIM,MISS_VAL,
**   +OUTPUT_BUFFER,ODIM,LENGTH)
**
** When calling PACKGRIB:
**   GRIB_LBL    is the user-filled array of codes and values to be packed into
**               the description sections of the grid
**   PDS_EXT     is the optional user-filled array of supplemental information
**               to be appended to the GRIB PDS
**   IPDIM       is the user-defined dimension of PDS_EXT
**               (use IPDIM = 0 for no inclusion of an extension to the GRIB
**               PDS)
**   GRIDPOINTS  is the user-filled array of real values that will be packed
**               into the grid.  GRIDPOINTS should be dimensioned at least NUM_I
**               by NUM_J, where NUM_I is the number of gridpoints in the
**               i-direction and NUM_J is the number of gridpoints in the
**               j-direction.
**   IGIDIM      the first dimension of GRIDPOINTS
**   MISS_VAL    is the value used to indicate missing data in the grid
**   ODIM        is the user-defined dimension of OUTPUT_BUFFER and should be at
**               least as large as the expected length of the GRIB grid
** 
** On return from PACKGRIB:
**   OUTPUT_BUFFER  is the buffer containing the GRIB representation of the grid
**   LENGTH         is the total length in octets (8-bit bytes) of the GRIB grid
** 
** Overview of GRIB_LBL:
**   The first 29 locations in GRIB_LBL contain information about the GRIB grid
**   that must always be present.  The remainder of GRIB_LBL contains the grid
**   description, which must also be present, but is dependent on the type of
**   grid being packed.
** 
** Description for GRIB_LBL(1) to GRIB_LBL(29):
**   GRIB_LBL(1):   GRIB Edition number
**   GRIB_LBL(2):   Length in octets (8-bit bytes) of the Product Definition
**                  Section
**   GRIB_LBL(3):   Parameter Table Version number
**   GRIB_LBL(4):   Center ID
**   GRIB_LBL(5):   Generating Process ID number
**   GRIB_LBL(6):   Grid Identification
**   GRIB_LBL(7):   0=Grid Description Section not included; 1=GDS included
**   GRIB_LBL(8):   Parameter Code
**   GRIB_LBL(9):   Level Type Code
**   GRIB_LBL(10):  Value for First Level
**   GRIB_LBL(11):  Value for Second Level (or 0)
**   GRIB_LBL(12):  Year (4-digits - YYYY)
**   GRIB_LBL(13):  Month
**   GRIB_LBL(14):  Day
**   GRIB_LBL(15):  Time (HHMM - HH=hour, MM=minutes)
**   GRIB_LBL(16):  Forecast Time Unit
**   GRIB_LBL(17):  P1 (or 0)
**   GRIB_LBL(18):  P2 (or 0)
**   GRIB_LBL(19):  Time Range Indicator
**   GRIB_LBL(20):  Number included in Average (or 0)
**   GRIB_LBL(21):  Number of Missing Grids in Average
**   GRIB_LBL(22):  Sub-center ID
**   GRIB_LBL(23):  Decimal scale factor
**   GRIB_LBL(24):  Binary Data Section flag
**   GRIB_LBL(25):  Width in bits of a packed data point
**   GRIB_LBL(26):  Length in octets of the Grid Description Section
**   GRIB_LBL(27):  Reserved - currently ignored
**   GRIB_LBL(28):  Reserved - currently ignored
**   GRIB_LBL(29):  Data representation type
** 
**   For Latitude/Longitude and Gaussian Lat/Lon Grids:
**     GRIB_LBL(30):  Number of points along a latitude circle
**     GRIB_LBL(31):  Number of points along a longitude meridian
**     GRIB_LBL(32):  Latitude of the first gridpoint (*1000)
**     GRIB_LBL(33):  Longitude of the first gridpoint (*1000)
**     GRIB_LBL(34):  Resolution and component flags
**     GRIB_LBL(35):  Latitude of the last gridpoint (*1000)
**     GRIB_LBL(36):  Longitude of the last gridpoint (*1000)
**     GRIB_LBL(37):  Latitude increment (*1000) for Lat/Lon grid
**                    -OR-
**                    Number of latitude circles between equator and pole for
**                    Gaussian Lat/Lon grid
**     GRIB_LBL(38):  Longitude increment (*1000)
**     GRIB_LBL(39):  Scanning mode flags
**
**   For Polar Stereographic Grids:
**     GRIB_LBL(30):  Number of points in the X-direction
**     GRIB_LBL(31):  Number of points in the Y-direction
**     GRIB_LBL(32):  Latitude of the first gridpoint (*1000)
**     GRIB_LBL(33):  Longitude of the first gridpoint (*1000)
**     GRIB_LBL(34):  Resolution and component flags
**     GRIB_LBL(35):  Longitude of grid orientation (*1000)
**     GRIB_LBL(36):  X-direction grid length in meters
**     GRIB_LBL(37):  Y-direction grid length in meters
**     GRIB_LBL(38):  Projection center flag
**     GRIB_LBL(39):  Scanning mode flags
** 
** 
** Description for PDS_EXT:
**   This array is free-form and can contain any values that the user wants to
**   specify, provided they are values that will fit into an 8-bit field.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void setBits(size_t *buf,size_t loc,size_t off,size_t bits)
{
  size_t mask=0;
  size_t size=sizeof(size_t)*8,wskip,bskip,lclear,rclear,left,right,more;
  size_t n;

/* no work to do */
  if (bits == 0)
    return;

/* create a mask to use when right-shifting (necessary because different
   compilers do different things when right-shifting a signed bit-field) */
  mask=1;
  for (n=1; n < size; n++) {
    mask<<=1;
    mask++;
  }

  if (bits > size) {
    fprintf(stderr,"Error: packing %d bits into a %d-bit field\n",bits,size);
    exit(1);
  }
  else {
/* get number of words and bits to skip before packing begins */
    wskip=off/size;
    bskip=off % size;
    lclear=bskip+bits;
    rclear=size-bskip;
    left= (rclear != size) ? (buf[wskip]&(mask<<rclear)) : 0;
    if (lclear <= size) {
/* all bits to be packed are in the current word */
/* clear the field to be packed */
	right=
        (lclear != size) ? (buf[wskip]&~(mask<<(size-lclear))) : 0;
/* fill the field to be packed */
	buf[wskip]=left|right|(loc<<(rclear-bits));
    }
    else {
/* bits to be packed cross a word boundary(ies) */
/* clear the bit field to be packed */
	more=bits-rclear;
	buf[wskip]=left|((loc>>more)&~(mask<<rclear));
/* clear the next (or part of the next) word and pack those bits */
	while (more > size) {
	  more-=size;
	  buf[++wskip]=((loc>>more)&~(mask<<size));
	}
	wskip++;
	more=size-more;
	right= (more != size) ? (buf[wskip]&~(mask<<more)) : 0;
	buf[wskip]=right|(loc<<more);
    }
  }
}

int round(double val)
{
  return ( (val >= 0.) ? (int)(val+0.5000001) : (int)(val-0.5000001) );
}

size_t real2ibm(double native_real)
{
  size_t ibm_real;
  size_t sign=0,fr=0;
  int exp=64;
  const double full=0xffffff;
  size_t size=sizeof(size_t)*8,off=0;

  if (native_real != 0.) {
    if (native_real < 0.) {
      sign=1;
      native_real=-native_real;
    }
    native_real/=pow(2.,-24.);
    while (exp > 0 && native_real < full) {
      native_real*=16.;
      exp--;
    }
    while (native_real > full) {
      native_real/=16.;
      exp++;
    }
    fr=round(native_real);

    if (size > 32) {
	off=size-32;
	setBits(&ibm_real,0,0,off);
    }
    setBits(&ibm_real,sign,off,1);
    setBits(&ibm_real,exp,off+1,7);
    setBits(&ibm_real,fr,off+8,24);
  }

  return ibm_real;
}

void packIS(int *grib_lbl,size_t *out_buf,size_t *out_buf_len,
  size_t *grib_length)
{
  size_t length;
  unsigned char *c_buf=(unsigned char *)out_buf;

  switch (grib_lbl[0]) {
    case 0:
	*grib_length=4;
	break;
    case 1:
	*grib_length=8;
/* GRIB Edition - set to 1 */
	c_buf[7]=1;
	break;
  }
  if (*grib_length > *out_buf_len) {
    fprintf(stderr,"Error: dimension of output buffer not large enough\n");
    exit(1);
  }
/* "GRIB" message */
  setBits(out_buf,0x47524942,0,32);
}

void packPDS(int *grib_lbl,size_t *pds_ext,size_t *pds_ext_length,
  size_t *out_buf,size_t *out_buf_len,size_t *grib_length,size_t *off)
{
  short sign;
  unsigned char flag=0x0;
  int dval,yr,cent,hr,min;
  size_t n;
  unsigned char *c_buf=(unsigned char *)out_buf;

  switch (grib_lbl[0]) {
    case 0:
	*off=32;
	*grib_length+=24;
/* length of PDS */
	setBits(out_buf,24,*off,24);
/* Edition number */
	setBits(out_buf,0,*off+24,8);
/* force the decimal scale factor to zero */
	grib_lbl[22]=0;
	break;
    case 1:
	*off=64;
	*grib_length+=grib_lbl[1];
/* length of PDS */
	setBits(out_buf,grib_lbl[1],*off,24);
/* table version */
	setBits(out_buf,grib_lbl[2],*off+24,8);
	break;
    default:
	fprintf(stderr,"Error: invalid GRIB edition number %d\n",grib_lbl[0]);
	exit(1);
  }

  if ( (grib_lbl[1] % 2) != 0) {
    fprintf(stderr,"Error: PDS length must be an even number of octets\n");
    exit(1);
  }

  if (*grib_length > *out_buf_len) {
    fprintf(stderr,"Error: dimension of output buffer not large enough\n");
    exit(1);
  }
/* center ID */
  setBits(out_buf,grib_lbl[3],*off+32,8);
/* generating process */
  setBits(out_buf,grib_lbl[4],*off+40,8);
/* grid type */
  setBits(out_buf,grib_lbl[5],*off+48,8);
/* flag to include GDS */
  if (grib_lbl[6] == 1) flag|=0x80;
  setBits(out_buf,flag,*off+56,8);
/* parameter code */
  setBits(out_buf,grib_lbl[7],*off+64,8);
/* level type */
  setBits(out_buf,grib_lbl[8],*off+72,8);
  switch (grib_lbl[8]) {
    case 100:
    case 103:
    case 105:
    case 107:
    case 109:
    case 111:
    case 113:
    case 115:
    case 125:
    case 160:
    case 200:
    case 201:
/* first level */
	setBits(out_buf,grib_lbl[9],*off+80,16);
	break;
    default:
/* first level */
	setBits(out_buf,grib_lbl[9],*off+80,8);
/* second level */
	setBits(out_buf,grib_lbl[10],*off+88,8);
  }
/* year */
  yr=grib_lbl[11] % 100;
  if (grib_lbl[0] > 0)
    if (yr == 0) yr=100;
  setBits(out_buf,yr,*off+96,8);
/* month */
  setBits(out_buf,grib_lbl[12],*off+104,8);
/* day */
  setBits(out_buf,grib_lbl[13],*off+112,8);
/* hour */
  hr=grib_lbl[14]*0.01;
  setBits(out_buf,hr,*off+120,8);
/* minutes */
  min=grib_lbl[14] % 100;
  setBits(out_buf,min,*off+128,8);
/* forecast units */
  setBits(out_buf,grib_lbl[15],*off+136,8);
/* P1 */
  setBits(out_buf,grib_lbl[16],*off+144,8);
/* P2 */
  setBits(out_buf,grib_lbl[17],*off+152,8);
/* time range */
  setBits(out_buf,grib_lbl[18],*off+160,8);
  switch (grib_lbl[18]) {
    case 3:
    case 4:
    case 51:
    case 113:
    case 114:
    case 115:
    case 116:
    case 117:
    case 123:
    case 124:
/* number in average */
      setBits(out_buf,grib_lbl[19],*off+168,16);
	break;
    default:
/* number in average */
      setBits(out_buf,0,*off+168,16);
  }
/* missing grids in average */
  setBits(out_buf,grib_lbl[20],*off+184,8);

/* if GRIB0, no more packing to do */
  if (grib_lbl[0] == 0) {
    *off+=192;
    return;
  }

/* century */
  cent=grib_lbl[11]*0.01;
  if (yr != 100) cent++;
  setBits(out_buf,cent,*off+192,8);
/* sub-center ID */
  setBits(out_buf,grib_lbl[21],*off+200,8);
/* decimal scale factor */
  dval=grib_lbl[22];
  sign=0;
  if (dval < 0) {
    sign=1;
    dval=-dval;
  }
  setBits(out_buf,sign,*off+208,1);
  setBits(out_buf,dval,*off+209,15);

  if (*pds_ext_length > 0) {
    if (grib_lbl[1] != 40+*pds_ext_length) {
	fprintf(stderr,"Error: bad PDS length %d\n",grib_lbl[1]);
	exit(1);
    }
    if ( (*pds_ext_length % 2) != 0) {
	fprintf(stderr,"Error: PDS supplement length must be an even number of "
        "octets\n");
	exit(1);
    }
    for (n=0; n < *pds_ext_length; n++)
	c_buf[48+n]=pds_ext[n];
  }
  *off+=grib_lbl[1]*8;
}

void packGDS(int *grib_lbl,size_t *out_buf,size_t *out_buf_len,
  size_t *grib_length,size_t *off)
{
  short sign;
  int dum;

  *grib_length+=grib_lbl[25];
  if (*grib_length > *out_buf_len) {
    fprintf(stderr,"Error: dimension of output buffer not large enough\n");
    exit(1);
  }
/* length of the GDS */
  setBits(out_buf,grib_lbl[25],*off,24);
/* data representation type */
  setBits(out_buf,grib_lbl[28],*off+40,8);
  switch (grib_lbl[28]) {
/* Latitude/Longitude grid */
    case 0:
/* Gaussian Lat/Lon grid */
    case 4:
/* number of latitudes */
	setBits(out_buf,grib_lbl[29],*off+48,16);
/* number of longitudes */
	setBits(out_buf,grib_lbl[30],*off+64,16);
/* latitude of first gridpoint */
	if (grib_lbl[31] < 0) {
	  sign=1;
	  grib_lbl[31]=-grib_lbl[31];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+80,1);
	setBits(out_buf,grib_lbl[31],*off+81,23);
/* longitude of first gridpoint */
	if (grib_lbl[32] < 0) {
	  sign=1;
	  grib_lbl[32]=-grib_lbl[32];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+104,1);
	setBits(out_buf,grib_lbl[32],*off+105,23);
/* resolution and component flags */
	setBits(out_buf,grib_lbl[33],*off+128,8);
/* latitude of last gridpoint */
	if (grib_lbl[34] < 0) {
	  sign=1;
	  grib_lbl[34]=-grib_lbl[34];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+136,1);
	setBits(out_buf,grib_lbl[34],*off+137,23);
/* longitude of last gridpoint */
	if (grib_lbl[35] < 0) {
	  sign=1;
	  grib_lbl[35]=-grib_lbl[35];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+160,1);
	setBits(out_buf,grib_lbl[35],*off+161,23);
/* longitude increment */
	setBits(out_buf,grib_lbl[37],*off+184,16);
/* latitude increment */
	setBits(out_buf,grib_lbl[36],*off+200,16);
/* scanning mode flag */
	setBits(out_buf,grib_lbl[38],*off+216,8);
	*off+=256;
	break;
/* Polar Stereographic grid */
    case 5:
/* number of points in the x-direction */
	setBits(out_buf,grib_lbl[29],*off+48,16);
/* number of points in the y-direction */
	setBits(out_buf,grib_lbl[30],*off+64,16);
/* latitude of first gridpoint */
	if (grib_lbl[31] < 0) {
	  sign=1;
	  grib_lbl[31]=-grib_lbl[31];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+80,1);
	setBits(out_buf,grib_lbl[31],*off+81,23);
/* longitude of first gridpoint */
	if (grib_lbl[32] < 0) {
	  sign=1;
	  grib_lbl[32]=-grib_lbl[32];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+104,1);
	setBits(out_buf,grib_lbl[32],*off+105,23);
/* resolution and component flags */
	setBits(out_buf,grib_lbl[33],*off+128,8);
/* longitude of grid orientation */
	if (grib_lbl[34] < 0) {
	  sign=1;
	  grib_lbl[34]=-grib_lbl[34];
	}
	else
	  sign=0;
	setBits(out_buf,sign,*off+136,1);
	setBits(out_buf,grib_lbl[34],*off+137,23);
/* x-direction grid length */
	setBits(out_buf,grib_lbl[35],*off+160,24);
/* y-direction grid length */
	setBits(out_buf,grib_lbl[36],*off+184,24);
/* projection center flag */
	setBits(out_buf,grib_lbl[37],*off+208,8);
/* scanning mode flag */
	setBits(out_buf,grib_lbl[38],*off+216,8);
	*off+=256;
	break;
  }
}

void packBDS(int *grib_lbl,float *gridpoints,size_t *gi_len,float miss_val,
  size_t *out_buf,size_t *out_buf_len,size_t *grib_length,size_t off)
{
  size_t n,m,*packed,ref_val;
  size_t bds_length,max_pack,num_missing=0;
  size_t num_points=grib_lbl[29]*grib_lbl[30],num_octets;
  int sign,eval,E=0,cnt=0,gcnt=0,bcnt=0;
  float max,min,range,e,d=pow(10.,grib_lbl[22]);
  size_t *bitmap,bms_length,missing_octets,ub;

  if ((grib_lbl[23] & 0x40) == 0) {
/* simple packing */
    num_octets=num_points*grib_lbl[24];
    switch (grib_lbl[5]) {
	case 21:
	case 22:
	case 23:
	case 24:
	case 25:
	case 26:
	case 61:
	case 62:
	case 63:
	case 64:
	  num_octets+=grib_lbl[24];
	  break;
    }
  }
  else {
/* second-order packing */
    fprintf(stderr,"Error: complex packing not currently supported\n");
    exit(1);
  }

  bitmap=(size_t *)malloc(sizeof(size_t)*num_points);
  cnt=0;
  min=1.e30;
  max=-min;
  for (n=0; n < grib_lbl[30]; n++) {
    for (m=0; m < *gi_len; m++) {
	if (m < grib_lbl[29]) {
	  if (gridpoints[gcnt] != miss_val) {
	    bitmap[cnt]=1;
	    if (gridpoints[gcnt] > max) max=gridpoints[gcnt];
	    if (gridpoints[gcnt] < min) min=gridpoints[gcnt];
	  }
	  else {
	    bitmap[cnt]=0;
	    num_missing++;
	  }
	  cnt++;
	}
	gcnt++;
    }
  }

  if (num_missing > 0) {
    bms_length=6+(num_points+7)/8;
    if ( (bms_length % 2) != 0) bms_length++;
    missing_octets=(grib_lbl[24]*num_missing+7)/8;
    if (missing_octets > bms_length) {
	num_octets-=(grib_lbl[24]*num_missing);

/* pack a bit-map section */
	*grib_length+=bms_length;
	if (*grib_length > *out_buf_len) {
	  fprintf(stderr,"Error: dimension of output buffer not large enough\n");
	  exit(1);
	}
/* set flag to show BMS included */
	out_buf[15]^=0x40;
/* length of the BMS */
	setBits(out_buf,bms_length,off,24);
/* unused bits at end of section */
	ub=((bms_length-6)*8) % num_points;
	setBits(out_buf,ub,off+24,8);
/* table reference set to show that a bitmap follows */
	setBits(out_buf,0,off+32,16);
	off+=48;
	for (n=0; n < num_points; n++) {
	  setBits(out_buf,bitmap[n],off,1);
	  off++;
	}
	off+=ub;
    }
    else {
/* let missing points become part of the scaling */
	if (miss_val > max) max=miss_val;
	if (miss_val < min) min=miss_val;
    }
  }

  if (min > 1.e29)
    min=0.;
  else {
    range=(max-min)*d;
    if (range == 0) {
	grib_lbl[24]=0;
	num_octets=0;
    }
    else {
	max_pack=pow(2.,grib_lbl[24])-1;
	if (range != 0.) {
	  while (round(range) <= max_pack) {
	    range*=2.;
	    E--;
	  }
	  while (round(range) > max_pack) {
	    range*=0.5;
	    E++;
	  }
	}
    }
  }

  num_octets=(num_octets+7)/8;
  bds_length=11+num_octets;
  if ((bds_length % 2) != 0) bds_length++;
  *grib_length+=bds_length;
  if (*grib_length > *out_buf_len) {
    fprintf(stderr,"Error: dimension of output buffer not large enough\n");
    exit(1);
  }
/* length of the BDS */
  setBits(out_buf,bds_length,off,24);
/* flag */
  setBits(out_buf,grib_lbl[23],off+24,4);
/* set the unused bits at the end of the BDS */
  if (grib_lbl[24] > 0)
    setBits(out_buf,(((bds_length-11)*8) % grib_lbl[24]),off+28,4);
  else
    setBits(out_buf,8,off+28,4);
/* binary scale factor */
  setBits(out_buf,abs(E),off+32,16);
  if (E < 0)
    setBits(out_buf,1,off+32,1);
  ref_val=real2ibm(min*d);
/* IBM representation of the reference value */
  setBits(out_buf,ref_val,off+48,32);
/* bit width of the packed data points */
  setBits(out_buf,grib_lbl[24],off+80,8);
  if (grib_lbl[24] == 0) {
    free(bitmap);
    return;
  }

  e=pow(2.,E);

  if ((grib_lbl[23] & 0x40) == 0) {
/* simple packing */
    off+=88;
    packed=(size_t *)malloc(sizeof(size_t)*num_points);
    cnt=0;
    switch (grib_lbl[28]) {
	case 0:
	case 4:
	  switch (grib_lbl[5]) {
	    case 23:
	    case 24:
	    case 26:
	    case 63:
	    case 64:
		setBits(out_buf,round((gridpoints[0]-min)*d/e),off,grib_lbl[24]);
		off+=grib_lbl[24];
		break;
	  }
	  gcnt=0;
	  for (n=0; n < grib_lbl[30]; n++) {
	    for (m=0; m < *gi_len; m++) {
		if (m < grib_lbl[29]) {
		  if (bitmap[bcnt] == 1) {
		    packed[cnt]=round((gridpoints[gcnt]-min)*d/e);
		    cnt++;
		  }
		  bcnt++;
		}
		gcnt++;
	    }
	  }
	  for (n=0; n < cnt; n++) {
	    setBits(out_buf,packed[n],off,grib_lbl[24]);
	    off+=grib_lbl[24];
	  }
	  switch (grib_lbl[5]) {
	    case 21:
	    case 22:
	    case 25:
	    case 61:
	    case 62:
		setBits(out_buf,
              round((gridpoints[*gi_len*grib_lbl[20]]-min)*d/e),off,
              grib_lbl[24]);
		break;
	  }
	  break;
	default:
	  fprintf(stderr,"Error: unable to pack grids with representation %d\n",
          grib_lbl[28]);
	  exit(1);
    }
  }

  free(bitmap);
  free(packed);
}

void packEND(size_t *out_buf,size_t *out_buf_len,size_t *grib_length)
{
  *grib_length+=4;
  if (*grib_length > *out_buf_len) {
    fprintf(stderr,"Error: dimension of output buffer not large enough\n");
    exit(1);
  }
/* "7777" */
  setBits(out_buf,0x37373737,(*grib_length-4)*8,32);
}

void packgrib_(int *grib_lbl,size_t *pds_ext,size_t *pds_ext_length,
  float *gridpoints,size_t *gi_len,float *miss_val,size_t *out_buf,
  size_t *out_buf_len,size_t *grib_length)
{
  size_t off;

  packIS(grib_lbl,out_buf,out_buf_len,grib_length);
  packPDS(grib_lbl,pds_ext,pds_ext_length,out_buf,out_buf_len,grib_length,&off);
  if (grib_lbl[6] == 1) packGDS(grib_lbl,out_buf,out_buf_len,grib_length,&off);
  packBDS(grib_lbl,gridpoints,gi_len,*miss_val,out_buf,out_buf_len,grib_length,
    off);
  packEND(out_buf,out_buf_len,grib_length);
  if (grib_lbl[0] == 1)
    setBits((size_t *)out_buf,*grib_length,32,24);
}