/*
** File: UNPACKGRIB.c
**
** Author:  Bob Dattore
**          NCAR/DSS
**          dattore@ucar.edu
**          (303) 497-1825
**
** Latest Revision:   2 Mar 99
**
** Purpose: to provide a single routine for unpacking 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 UNPACKGRIB.
**
**          2) The user is expected to have some understanding of the GRIB
**             format, as some of the information returned by UNPACKGRIB
**             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) UNPACKGRIB has been tested on a Cray running UNICOS - 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  INPUT_BUFFER(ODIM)
** 
**    CALL UNPACKGRIB(GRIB_LBL,PDS_EXT,IPDIM,GRIDPOINTS,IGIDIM,MISS_VAL,
**   +INPUT_BUFFER,IDIM,LENGTH)
** 
** When calling UNPACKGRIB:
**   IPDIM         is the user-defined dimension of PDS_EXT
**   IGIDIM        the first dimension of GRIDPOINTS
**                 to unpack only header information (GRIB_LBL(1-39)), set
**                 IGIDIM to 0
**   MISS_VAL      is the value to use to indicate missing data in the grid
**   INPUT_BUFFER  is the buffer containing the GRIB representation of the grid
**   IDIM          is the user-defined dimension of INPUT_BUFFER
**
** On return from UNPACKGRIB:
**   GRIB_LBL    is the array of codes and values that have been unpacked from
**               the description sections of the grid
**   PDS_EXT     is the array of supplemental information, if it exists, that
**               was found to be appended to the GRIB PDS
**   IPDIM       = 0 if no extension was found; otherwise is the length in
**               octets of the extension that was found
**   GRIDPOINTS  is the array of real values that has been unpacked from the
**               grid.  GRIDPOINTS should be dimensioned at least as large as
**               the expected grid size.
**   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.
**   The remainder of GRIB_LBL contains the grid description, which is
**   dependent on the type of grid that was unpacked.
** 
** 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;
**                  -1=Grid type is not recognized (GRIB_LBL(30) through
**                  GRIB_LBL(39) should be ignored, and GRIDPOINTS will contain
**                  the gridpoint data as a single stream)
**   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 contains any 8-bit values that were found after
**   the end of the standard PDS, but before the next GRIB seciton.
*/

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

void getBits(size_t *buf,int *loc,size_t off,size_t bits)
{
  size_t mask=0;
  size_t size=sizeof(size_t)*8,wskip;
  int rshift;
  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: unpacking %d bits into a %d-bit field\n",bits,size);
    exit(1);
  }
  else {
/* get number of words to skip before unpacking begins */
    wskip=off/size;
/* right shift the bits in the packed buffer word to eliminate unneeded bits */
    rshift=size-(off % size)-bits;
/* check for a packed field spanning two words */
    if (rshift < 0) {
	*loc=(buf[wskip]<<-rshift);
	*loc+=(buf[++wskip]>>(size+rshift))&~(mask<<-rshift);
    }
    else
	*loc=(buf[wskip]>>rshift);
/* remove any unneeded leading bits */
    if (bits != size) *loc&=~(mask<<bits);
  }
}

double ibm2real(size_t *buf,size_t off)
{
  int sign,exp,fr;
  double native_real;

  getBits(buf,&sign,off,1);
  getBits(buf,&exp,off+1,7);
  exp-=64;
  getBits(buf,&fr,off+8,24);
  native_real=pow(2.,-24.)*(double)fr*pow(16.,(double)exp);

  return (sign == 1) ? -native_real : native_real;
}

void unpackIS(int *grib_lbl,size_t *in_buf,size_t *in_buf_len,
  size_t *grib_length)
{
  char message[5];

  memcpy(message,(char *)in_buf,4);  /* "GRIB" message */
  message[4]='\0';
  if (strcmp(message,"GRIB") != 0) {
    fprintf(stderr,"Error: not a GRIB record\n");
    exit(1);
  }
  getBits(in_buf,(int *)grib_length,32,24);
  if (*grib_length == 24) {
    grib_lbl[0]=0;
    grib_lbl[1]=*grib_length;

/* add the four bytes for 'GRIB' + 3 bytes for the length of the section
** following the PDS */
    *grib_length+=7;
  }
  else
    grib_lbl[0]=1;

  if (*grib_length > *in_buf_len) {
    fprintf(stderr,"Error: user buffer too small for GRIB length %d\n",
      *grib_length);
    exit(1);
  }
}

void unpackPDS(int *grib_lbl,int *pds_ext,size_t *pds_ext_length,size_t *in_buf,
  size_t *in_buf_len,size_t *off,size_t *grib_length)
{
  short ext_length;
  size_t n;
  int flag,min,cent,sign;
  char *c_buf=(char *)in_buf;

  if (grib_lbl[0] == 0)
    *off=32;
  else {
    *off=64;
    getBits(in_buf,&grib_lbl[1],*off,24); /* length of PDS */
    getBits(in_buf,&grib_lbl[2],*off+24,8);  /* table version */
  }
  getBits(in_buf,&grib_lbl[3],*off+32,8);  /* center ID */
  getBits(in_buf,&grib_lbl[4],*off+40,8);  /* generating process */
  getBits(in_buf,&grib_lbl[5],*off+48,8);  /* grid type */
  getBits(in_buf,&flag,*off+56,8);
  if ( (flag & 0x80) == 0x80)  /* indication of GDS */
    grib_lbl[6]=1;
  else
    grib_lbl[6]=0;
  if ( (flag & 0x40) == 0x40)  /* indication of BMS */
    grib_lbl[6]+=10;
  getBits(in_buf,&grib_lbl[7],*off+64,8);  /* parameter code */
  getBits(in_buf,&grib_lbl[8],*off+72,8);  /* level type */
  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:
	getBits(in_buf,&grib_lbl[9],*off+80,16);  /* first level */
	grib_lbl[10]=0;
	break;
    default:
	getBits(in_buf,&grib_lbl[9],*off+80,8);  /* first level */
	getBits(in_buf,&grib_lbl[10],*off+88,8);  /* second level */
  }
  getBits(in_buf,&grib_lbl[11],*off+96,8);  /* year of the century */
  getBits(in_buf,&grib_lbl[12],*off+104,8);  /* month */
  getBits(in_buf,&grib_lbl[13],*off+112,8);  /* day */
  getBits(in_buf,&grib_lbl[14],*off+120,8);  /* hour */
  getBits(in_buf,&min,*off+128,8);  /* minutes */
  grib_lbl[14]=grib_lbl[14]*100+min;
  getBits(in_buf,&grib_lbl[15],*off+136,8);  /* forecast units */
  getBits(in_buf,&grib_lbl[16],*off+144,8);  /* P1 */
  getBits(in_buf,&grib_lbl[17],*off+152,8);  /* P2 */
  getBits(in_buf,&grib_lbl[18],*off+160,8);  /* time range */
  switch (grib_lbl[18]) {
    case 3:
    case 4:
    case 51:
    case 113:
    case 114:
    case 115:
    case 116:
    case 117:
    case 123:
    case 124:
      getBits(in_buf,&grib_lbl[19],*off+168,16);  /* number in average */
	break;
    default:
	grib_lbl[19]=0;  /* number in average */
  }
  getBits(in_buf,&grib_lbl[20],*off+184,8);  /* missing grids in average */

/* if GRIB0, done with PDS */
  if (grib_lbl[0] == 0) {
    *pds_ext_length=0;
    *off+=192;
    return;
  }

  getBits(in_buf,&cent,*off+192,8);  /* century */
  grib_lbl[11]+=(cent-1)*100;
  getBits(in_buf,&grib_lbl[21],*off+200,8);  /* sub-center ID */
  getBits(in_buf,&sign,*off+208,1);
  getBits(in_buf,&grib_lbl[22],*off+209,15);  /* decimal scale factor */
  if (sign == 1)
    grib_lbl[22]=-grib_lbl[22];

  *off+=224;
  if (grib_lbl[1] > 28) {
    if (grib_lbl[1] < 40) {
	fprintf(stderr,"Warning: PDS extension is in wrong location\n");
	ext_length=grib_lbl[1]-28;
	if (ext_length > *pds_ext_length) {
	  fprintf(stderr,"Error: buffer not large enough to hold the PDS "
          "extension\n");
	  exit(1);
	}
	*pds_ext_length=ext_length;
	for (n=0; n < *pds_ext_length; n++)
	  pds_ext[n]=c_buf[36+n];
	*off+=*pds_ext_length*8;
    }
    else {
	ext_length=grib_lbl[1]-40;
	if (ext_length > *pds_ext_length) {
	  fprintf(stderr,"Error: buffer not large enough to hold the PDS "
          "extension\n");
	  exit(1);
	}
	*pds_ext_length=ext_length;
	for (n=0; n < *pds_ext_length; n++)
	  pds_ext[n]=c_buf[48+n];
	*off+=(*pds_ext_length+12)*8;
    }
  }
  else
    *pds_ext_length=0;
}

void unpackGDS(int *grib_lbl,size_t *in_buf,size_t *in_buf_len,size_t *off,
  size_t *grib_length)
{
  int sign,dum;

  getBits(in_buf,&grib_lbl[25],*off,24);  /* length of the GDS */

/* test the length of the user buffer for GRIB0 */
  if (grib_lbl[0] == 0) {
    *grib_length+=grib_lbl[25];
    if (*grib_length > *in_buf_len) {
	fprintf(stderr,"Error: user buffer too small for GRIB length %d\n",
        *grib_length);
	exit(1);
    }
  }
  getBits(in_buf,&grib_lbl[28],*off+40,8);  /* data representation type */
  switch (grib_lbl[28]) {
/* Latitude/Longitude grid */
    case 0:
/* Gaussian Lat/Lon grid */
/* Rotated Lat/Lon grid */
    case 4:
    case 10:
	getBits(in_buf,&grib_lbl[29],*off+48,16);  /* number of latitudes */
	getBits(in_buf,&grib_lbl[30],*off+64,16);  /* number of longitudes */
	getBits(in_buf,&sign,*off+80,1);
	getBits(in_buf,&grib_lbl[31],*off+81,23);  /* latitude of first
                                                    gridpoint */
	if (sign == 1)
	  grib_lbl[31]=-grib_lbl[31];
	getBits(in_buf,&sign,*off+104,1);
	getBits(in_buf,&grib_lbl[32],*off+105,23);  /* longitude of first
                                                     gridpoint */
	if (sign == 1)
	  grib_lbl[32]=-grib_lbl[32];
	getBits(in_buf,&grib_lbl[33],*off+128,8);  /* resolution and component
                                                    flags */
	getBits(in_buf,&sign,*off+136,1);
	getBits(in_buf,&grib_lbl[34],*off+137,23);  /* latitude of last
                                                     gridpoint */
	if (sign == 1)
	  grib_lbl[34]=-grib_lbl[34];
	getBits(in_buf,&sign,*off+160,1);
	getBits(in_buf,&grib_lbl[35],*off+161,23);  /* longitude of last
                                                     gridpoint */
	if (sign == 1)
	  grib_lbl[35]=-grib_lbl[35];
	getBits(in_buf,&grib_lbl[37],*off+184,16);  /* longitude increment */
	getBits(in_buf,&grib_lbl[36],*off+200,16);  /* latitude increment */
	getBits(in_buf,&grib_lbl[38],*off+216,8);  /* scanning mode flag */
	*off+=256;
	break;
/* Polar Stereographic grid */
    case 5:
	getBits(in_buf,&grib_lbl[29],*off+48,16);  /* number of x-points */
	getBits(in_buf,&grib_lbl[30],*off+64,16);  /* number of y-points */
	getBits(in_buf,&sign,*off+80,1);
	getBits(in_buf,&grib_lbl[31],*off+81,23);  /* latitude of first
                                                    gridpoint */
	if (sign == 1)
	  grib_lbl[31]=-grib_lbl[31];
	getBits(in_buf,&sign,*off+104,1);
	getBits(in_buf,&grib_lbl[32],*off+105,23);  /* longitude of first
                                                     gridpoint */
	if (sign == 1)
	  grib_lbl[32]=-grib_lbl[32];
	getBits(in_buf,&grib_lbl[33],*off+128,8);  /* resolution and component
                                                    flags */
	getBits(in_buf,&sign,*off+136,1);
	getBits(in_buf,&grib_lbl[34],*off+137,23);  /* longitude of grid
                                                     orientation */
	if (sign == 1)
	  grib_lbl[34]=-grib_lbl[34];
	getBits(in_buf,&grib_lbl[35],*off+160,24);  /* x-direction grid length */
	getBits(in_buf,&grib_lbl[36],*off+184,24);  /* y-direction grid length */
	getBits(in_buf,&grib_lbl[37],*off+208,8);  /* projection center flag */
	getBits(in_buf,&grib_lbl[38],*off+216,8);  /* scanning mode flag */
	*off+=256;
	break;
  }
}

void unpackBDS(int *grib_lbl,float *gridpoints,size_t *gi_len,float miss_val,
  size_t *in_buf,size_t *in_buf_len,size_t off,size_t *grib_length)
{
  size_t n,m;
  size_t boff,bcnt=0;
  size_t num_packed;
  int bms_length,bds_length,sign,ub;
  int tref,*bitmap=NULL;
  int E,*packed,cnt,gcnt=0;
  float ref_val,e,d=pow(10.,grib_lbl[22]);

  if (grib_lbl[6] >= 10) {
    getBits(in_buf,&bms_length,off,24);

/* test the length of the user buffer for GRIB0 */
    if (grib_lbl[0] == 0) {
	*grib_length+=bms_length;
	if (*grib_length > *in_buf_len) {
	  fprintf(stderr,"Error: user buffer too small for GRIB length %d\n",
          *grib_length);
	  exit(1);
	}
    }
    getBits(in_buf,&ub,off+24,8);
    getBits(in_buf,&tref,off+32,16);
    if (tref != 0) {
	fprintf(stderr,"Error: unknown pre-defined bit-map %d\n",tref);
	exit(1);
    }
    num_packed=(bms_length-6)*8-ub;
    bitmap=(int *)malloc(sizeof(int)*num_packed);
    boff=off+48;
    for (n=0; n < num_packed; n++) {
	getBits(in_buf,&bitmap[n],boff,1);
	boff++;
    }
    off+=bms_length*8;
    grib_lbl[6]%=10;
  }

  getBits(in_buf,&bds_length,off,24);  /* length of the BDS */

/* test the length of the user buffer for GRIB0 */
  if (grib_lbl[0] == 0) {
    *grib_length+=(bds_length+1);
    if (*grib_length > *in_buf_len) {
	fprintf(stderr,"Error: user buffer too small for GRIB length %d\n",
        *grib_length);
	exit(1);
    }
  }
  getBits(in_buf,&grib_lbl[23],off+24,4);  /* flag */
  getBits(in_buf,&ub,off+28,4);
  getBits(in_buf,&grib_lbl[24],off+80,8);  /* bit width of the packed data
                                              points */

/* no unpacking of gridpoint data */
  if (*gi_len == 0)
    return;

  getBits(in_buf,&sign,off+32,1);
  getBits(in_buf,&E,off+33,15); /* binary scale factor */
  if (sign == 1)
    E=-E;
  e=pow(2.,E);
  ref_val=ibm2real(in_buf,off+48);  /* reference value */

  if ((grib_lbl[23] & 0x40) == 0) {
/* simple packing */
    off+=88;
    num_packed=(bds_length*8-88-ub)/grib_lbl[24];
    packed=(int *)malloc(sizeof(int)*num_packed);
    cnt=0;
    switch (grib_lbl[28]) {
/* Latitude/Longitude grid */
	case 0:
/* Gaussian Lat/Lon grid */
/* Rotated Lat/Lon grid */
	case 4:
	case 10:
	  switch (grib_lbl[5]) {
	    case 23:
	    case 24:
	    case 26:
	    case 63:
	    case 64:
		off+=grib_lbl[24];
		break;
	  }
/* Polar Stereographic grid */
	case 5:
	  for (n=0; n < num_packed; n++) {
	    getBits(in_buf,&packed[n],off,grib_lbl[24]);
	    off+=grib_lbl[24];
	  }

/* no recognized GDS, so just unpack the stream of gridpoints */
	  if (grib_lbl[29] == 0 && grib_lbl[30] == 0) {
	    for (n=0; n < num_packed; n++) {
		if (bitmap == NULL || (bitmap != NULL && bitmap[bcnt++] == 1))
		  gridpoints[n]=(ref_val+packed[n]*e)/d;
		else
		  gridpoints[n]=miss_val;
	    }
	  }
	  else {
	    gcnt=0;
	    for (n=0; n < grib_lbl[30]; n++) {
		for (m=0; m < *gi_len; m++) {
		  if (m < grib_lbl[29]) {
		    if (bitmap == NULL || (bitmap != NULL && bitmap[bcnt++] == 1))
			gridpoints[gcnt]=(ref_val+packed[cnt++]*e)/d;
		    else
			gridpoints[gcnt]=miss_val;
		  }
		  gcnt++;
		}
	    }
	  }
	  break;
	default:
	  fprintf(stderr,"Error: unable to unpack grids with representation %d\n",
          grib_lbl[28]);
	  exit(1);
    }
  }
  else {
/* second-order packing */
    fprintf(stderr,"Error: complex packing not currently supported\n");
    exit(1);
  }

  free(packed);
}

void unpackEND(size_t *in_buf,size_t *in_buf_len,size_t grib_length)
{
  char *c_buf=(char *)in_buf,message[5];

  memcpy(message,&c_buf[grib_length-4],4);  /* "7777" end message */
  message[4]='\0';
  if (strcmp(message,"7777") != 0)
    fprintf(stderr,"Warning: no end section found\n");
}

void UNPACKGRIB(int *grib_lbl,int *pds_ext,size_t *pds_ext_length,
  float *gridpoints,size_t *gi_len,float *miss_val,size_t *in_buf,
  size_t *in_buf_len,size_t *grib_length)
{
  size_t n,off;

  for (n=0; n < 39; n++)
    grib_lbl[n]=0;

  unpackIS(grib_lbl,in_buf,in_buf_len,grib_length);
  unpackPDS(grib_lbl,pds_ext,pds_ext_length,in_buf,in_buf_len,&off,grib_length);
  if ((grib_lbl[6] % 10) == 1) unpackGDS(grib_lbl,in_buf,in_buf_len,&off,
    grib_length);
  unpackBDS(grib_lbl,gridpoints,gi_len,*miss_val,in_buf,in_buf_len,off,
    grib_length);
  unpackEND(in_buf,in_buf_len,*grib_length);
}