/* This routine will unpack a GRIB 1 format file and put all of the grids into
   an AIX McIDAS type gridfile. No grid interpolation is done */
#include <stdio.h>
#include <string.h>
#include <fcntl.h>
#include <math.h>
#include <fcntl.h>
#include <unistd.h>
#define MAIN_R 1
#include "ktypes.h"
#include "GribHead.h"
#define GRID 6000

void etoa (char *, int);
long bword (char *, int , int );
long exbits (int , int , char *);
int  x9ie32 (float *input, float *output, int niter);
int savegd (int, int);
int lscale (unsigned short int);

float unpack[15000];
long docsec[1000];	/* grib record documentation blocks */
long title[9];		/* grid title */
char out_file[140];
int MaxGrd;		/* grids in gridfile */
int Min,Max;
int bug = 0;
int Scale;

void main (int argc, char **argv)
{
 long grid[GRID];
 long blank = 0x20202020;	/* blanks */
 int fdi, lpos, start_byte, n_bytes, n_read, i,j,ig;
 int n_grids;		/* number of grids, from header */
 const int lenis = 8;	/* length of IS Indicator Section */
 int lenpds;		/* length of Product Definition Section in bytes */
 int lengds, lenbms, lenbds;
 int locpds, locgds, locbms, locbds;
 int grid_id;		/* grid identification byte 7 */
			/* GDS Grid Descriptor File */
 int flag_gds_bms;	/* presence of GDS or BMS */
 int i_type,i_level;
 unsigned short int dscale;		/* decimal scale fcator */
 unsigned short int bscale;		/* binary scale fcator */
 int scaleb;		/* binary scaling fcator */
 int scaled;		/* decimal scaling factor */
 union {
  long  i;		/* reference value */
  float f;} ref_val;
 int n_bits;		/* number of bits */
 int exw, is;
 float fexw, fscb,fscd;	/* float for use in unpacking grid */
 float q9f1,q9f2;	/* for call to q9ie32 */
 int qstat,qno;		/* for call to q9ie32 */
 int end_flag = 0;	/* end of file flag for call to savegd */
 int kt;		/* for ktype (index) */
 int beggrid, endgrid;	/* beginning, ending grid to load into gridfile */
 char cfile[140];
 char cline[140];	/* read buffer */
 FILE *fpc;

 printf("gribx1 version 29 Sep 93   1\n");
 if (argc == 1) goto help;
 if (!strncmp(argv[1],"help",4)) {goto help;}
 for (i=0;i<8;i++) {title[i]=blank;} title[8] = 0;
 /* The first parameter in command line is the name of the command file */
 /* from which the run parameters are taken from.  The first parameter is */
 /* the name of the GRIB file (followed by a line feed), the second */
 /* parameter is the output gridfile name (lf), third is the maximum number */
 /* of grids to reserve for the gridfile if it is being created. */
 if (argc == 2)
  {
   printf("file name: %s \n",argv[1]) ;
   fpc = fopen(argv[1],"r");
   if (fgets(cfile,140,fpc) == NULL)
    {printf("input data set not given\n"); return;}
   i = strlen (cfile);
   cfile[i-1] = '\0';  /* blot out the line feed character */
   printf ("input file: %s \n",cfile);
   if (fgets(out_file,140,fpc) == NULL)
    {printf("output data set not given\n"); return;}
   i = strlen (out_file);
   out_file[i] = '\0'; out_file[i-1] = '\0';
   printf ("output file: %s \n",out_file);
   /* Get size of output file 3rd parameter*/
   if (fgets(cline,140,fpc) == NULL)
    {printf("number of grids not given\n"); MaxGrd = -1;}
    else
    {i = strlen (cline); cline[i-1] = '\0'; sscanf(cline,"%d",&MaxGrd);}
   printf("Max Grids in gridfile: %d\n",MaxGrd);
   /* Beginning and ending grid number 0-based, blank line- use defaults */
   if (fgets(cline,140,fpc) == NULL)
    {printf("default begin & end grid of 0 & number of grids in file\n");
     beggrid = 0; endgrid = 9999;}
    else
    {sscanf(cline,"%d %d",&beggrid,&endgrid);
     printf("begin & end grid: %d %d \n",beggrid,endgrid);}
   /* Title */
   if (fgets((char *)title,33,fpc) == NULL)
    {printf("no title\n");}
    else
    {for(i=0;i<32;i++) {if( *(((char *)&title)+i) == '\n')
                         {*(((char *)&title)+i) = ' '; }
                        if( *(((char *)&title)+i) == '\0')
                         {*(((char *)&title)+i) = ' '; } }
     printf("title: %s<\n",title);}
   fclose(fpc);
  }
  else
  {
   if (argc < 2) 
    {printf("need at least 2 arguements\n"); goto help;}
   strncpy(cfile,argv[1],139);
   printf ("input file: %s \n",cfile);
   strncpy(out_file,argv[2],139);
   printf ("output file: %s \n",out_file);
   if (argc >= 4)
    {sscanf(argv[3],"%d",&MaxGrd);}
    else
    {MaxGrd = 250;}
   printf("Max Grids in gridfile: %d\n",MaxGrd);
   if (argc >= 5)
    {sscanf(argv[4],"%d",&beggrid);}
    else
    {beggrid = 0;}
   printf("beginning grid %d \n",beggrid);
   if (argc >= 6)
    {sscanf(argv[5],"%d",&endgrid);}
    else
    {endgrid = 9999;}
   printf("ending grid %d \n",endgrid);
   if (argc >= 7)
    {strncpy((char *)title,argv[6],32);
     printf("title: %s<\n",title);}
  }
 /* Open the input file containing NMC grib data with read only access */
 fdi = open (cfile, O_RDONLY, 0) ;
 if (fdi < 0) {printf("open error\n") ; goto error ;}
 printf ("fdi= %8d \n",fdi);
 start_byte = 0;  n_bytes = 100;
 if (lpos = lseek(fdi, start_byte, 0) < 0)
    {printf ("lseek error\n") ; goto error ; }
 n_read = read (fdi, &header, n_bytes ) ;
 if (n_read != n_bytes) {printf("read error %d %d \n",n_read,n_bytes);}
 etoa ((char *)header.label,32);
 printf ("header record \n");
 printf ("label: %.32s \n",header.label);
 printf ("length of index: %5d \n",header.len_index);
 printf ("number of grids: %5d \n",header.no_grids);
 n_grids = header.no_grids;
 if (MaxGrd < 1) MaxGrd = n_grids;	/* reserve for maxgrids of gridfile */
 printf ("MaxGrd %d \n",MaxGrd);
 if (endgrid >= n_grids) endgrid = n_grids - 1;	/* reset if endgrid too big */
 /* read in the index */
 n_bytes = 2 * header.len_index;
 start_byte = 100;
 if (lpos = lseek(fdi, start_byte, 0) < 0)
    {printf ("lseek error\n") ; goto error ; }
 n_read = read (fdi, indx, n_bytes ) ;
 if (n_read != n_bytes) {printf("read error %d %d \n",n_read,n_bytes);}
 for (i=0; i<header.no_grids&&i<5; i++)
  {printf("%08x ",indx[i].loc);
   for (j=0; j<7; j++) {printf("%08x ",indx[i].pds[j]);}
   printf("%08x\n",indx[i].len);}
 printf ("\n");  /* final line skip before data records */

 /* read in a data record */
 for (ig=0; ig<n_grids; ig++)
  {start_byte = indx[ig].loc-1;	/* location byte is 1 based */
   n_bytes = indx[ig].len;
   if (ig < beggrid || ig > endgrid) continue;
   if (lpos = lseek(fdi, start_byte, 0) < 0)
      {printf ("lseek error\n") ; goto error ; }
   if (n_bytes > GRID*4)
      {printf("n bytes = %d GRID = %d start byte %d \n",n_bytes,GRID,start_byte);
       continue; }
   n_read = read (fdi, grid, n_bytes ) ;
   if (n_read != n_bytes) {printf("read error %d %d \n",n_read,n_bytes);}
   j = n_bytes/4;
   for (i=0; i<8; i++)
      printf("%08x%c", grid[i], (i%8 == (8-1) || i==j-1) ? '\n' : ' ');
   locpds = lenis;
   lenpds = bword ( (char *)grid, locpds, 3);
   grid_id = bword ( (char *)grid, locpds+6, 1);
   flag_gds_bms = bword ( (char *)grid, locpds+7, 1);
   i_type  = bword ( (char *)grid, locpds+8, 1);
   i_level = bword ( (char *)grid, locpds+9, 1);
   if (flag_gds_bms != 0) {printf("** GDS or BMS present %08x \n",flag_gds_bms);
 			continue; }
   dscale = bword ( (char *)grid, locpds+26, 2);
   if (lenpds > GRID*4) {printf("lenpds too large %d\n",lenpds); continue;}
   lengds = 0; lenbms = 0;	/* TEMPORARY */
   locgds = locpds + lenpds;
   if ((flag_gds_bms & 8) > 0)
    {lengds = bword ( (char *)grid, locgds, 3); printf("GDS present\n");}
    else
    {lengds = 0;}	/* grid descriptor section not present */
   locbms = locgds + lengds;
   locbds = locbms + lenbms;
   lenbds = bword ( (char *)grid, locbds, 3);
   for (i=0; i<(lenpds+lengds+lenbms); i++) {docsec[i]=grid[i];}
   bscale = bword ( (char *)grid, locbds+4, 2);
   ref_val.i = bword ( (char *)grid, locbds+6, 4);
   q9f1 = ref_val.f; qno = 1;
   /* q9ie32 (&q9f1,&q9f2,&qno,&qstat); */
   qstat = x9ie32 (&ref_val.f,&q9f2,qno);
   if (qstat != 0) {printf("q9ie32 error \n"); continue;}
   ref_val.f = q9f2;
   n_bits = bword ( (char *)grid, locbds+10, 1);
   scaled = lscale(dscale);
   scaleb = lscale(bscale);
   printf ("lenbds %7d n bits %2d bscale %3d dscale %hd ",
           lenbds,n_bits,scaleb,scaled);
   Scale = scaled;
   if (n_bits == 0) {printf("n_bits = 0 SKIP GRIB\n"); continue;}
   /* Check for recognized ktypes. Do not unpack unrecognized ktypes */
   for(kt=0; kt<PROJSIZ; kt++) {if(grid_id == iptype[kt]) goto GoodGrid;}
   printf ("ktype %d not recognized SKIP GRIB\n",grid_id);
   continue;	/* do not process this grid */
   GoodGrid:
   printf("unpack grid ktype %d size %d seq no %d\n",kt,ipgsiz[kt],ig);
   fflush (NULL);
   Min = 99999999;  Max = -Min;
   if(scaleb != 0) {fscb = pow(2.0,(double)scaleb);} else {fscb = 1.0;}
   if(scaled != 0) {fscd = pow(10.,(double)scaled);} else {fscd = 1.0;}
   is = (locbds + 11) * 8;
   for (i=0; i<ipgsiz[kt]; i++)
    {
     exw = exbits(is,n_bits,(char *)grid); is = is + n_bits;
     fexw = (float)exw;
     if (scaleb != 0) {fexw = fexw * fscb;}
     fexw = fexw + ref_val.f;
     exw = (int)fexw;	/* should do scaling here, if any */
     unpack[i] = exw;
     if (exw < Min) Min = exw;  if (exw > Max) Max = exw;
    }
   savegd (end_flag,grid_id);
   printf ("\n");
  }		/* end for on ig for grids */
 end_flag = 1;		/* signal end of file */
 savegd (end_flag,grid_id);
 return;

help:
 printf("gribx1: unpacks all gribfile grids and puts them into a gridfile\n");
 printf("There are 2 ways to enter parameters:\n");
 printf("command line: gribx1 (input file) (output gridfile) (max grids) "
        "(beginning grid 0-based) (ending grid) (\"title\")\n");
 printf("command file: gribx1 (command file name)\n");
 printf("format of command file:\n");
 printf(" line 1: name of the grib file - input (lf)\n");
 printf(" line 2: name of the gridfile - output (lf)\n");
 printf(" line 3: maximum number of grids in gridfile\n");
 printf(" line 4: starting grid (def 0) ending grid (def 9999)\n");
 printf(" line 5: title 8 words \n");
 return;
error: return;
}