/src/gdal/frmts/grib/degrib/g2clib/comunpack.c

Source
#include <float.h>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include "grib2.h"

#ifndef DoubleToFloatClamp_defined
#define DoubleToFloatClamp_defined
static float DoubleToFloatClamp(double val) {
   if (val >= FLT_MAX) return FLT_MAX;
   if (val <= -FLT_MAX) return -FLT_MAX;
   return (float)val;
}
#endif

int comunpack(unsigned char *cpack,g2int cpack_length,g2int lensec,g2int idrsnum,g2int *idrstmpl,g2int ndpts,g2float *fld)
////$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    comunpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-10-29
//
// ABSTRACT: This subroutine unpacks a data field that was packed using a
//   complex packing algorithm as defined in the GRIB2 documentation,
//   using info from the GRIB2 Data Representation Template 5.2 or 5.3.
//   Supports GRIB2 complex packing templates with or without
//   spatial differences (i.e. DRTs 5.2 and 5.3).
//
// PROGRAM HISTORY LOG:
// 2002-10-29  Gilbert
// 2004-12-16  Gilbert  -  Added test ( provided by Arthur Taylor/MDL )
//                         to verify that group widths and lengths are
//                         consistent with section length.
//
// USAGE:    int comunpack(unsigned char *cpack,g2int lensec,g2int idrsnum,
//                         g2int *idrstmpl, g2int ndpts,g2float *fld)
//   INPUT ARGUMENT LIST:
//     cpack    - pointer to the packed data field.
//     lensec   - length of section 7 (used for error checking).
//     idrsnum  - Data Representation Template number 5.N
//                Must equal 2 or 3.
//     idrstmpl - pointer to the array of values for Data Representation
//                Template 5.2 or 5.3
//     ndpts    - The number of data values to unpack
//
//   OUTPUT ARGUMENT LIST:
//     fld      - Contains the unpacked data values.  fld must be allocated
//                with at least ndpts*sizeof(g2float) bytes before
//                calling this routine.
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:
//
//$$$//
{

      g2int nbitsd=0,isign;
      g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
      g2int *ifld=NULL,*ifldmiss=NULL;
      g2int *gref=NULL,*gwidth=NULL,*glen=NULL;
      g2int itype,ngroups,nbitsgref,nbitsgwidth,nbitsglen;
      g2int msng1,msng2;
      g2float ref,bscale,dscale,rmiss1,rmiss2;
      g2int totBit, totLen;
      g2int errorOccurred = 0;

      //printf('IDRSTMPL: ',(idrstmpl(j),j=1,16)
      rdieee(idrstmpl+0,&ref,1);
//      printf("SAGTref: %f\n",ref);
      bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
      dscale = DoubleToFloatClamp(int_power(10.0,-idrstmpl[2]));
      nbitsgref = idrstmpl[3];
      itype = idrstmpl[4];
      ngroups = idrstmpl[9];
      nbitsgwidth = idrstmpl[11];
      nbitsglen = idrstmpl[15];
      if (idrsnum == 3)
         nbitsd=idrstmpl[17]*8;

      //   Constant field

      if (ngroups == 0) {
         for (j=0;j<ndpts;j++) fld[j]=ref;
         return(0);
      }

      /* To avoid excessive memory allocations. Not completely sure */
      /* if this test is appropriate (the 10 and 2 are arbitrary), */
      /* but it doesn't seem to make sense to have ngroups much larger than */
      /* ndpts */
      if( ngroups < 0 || ngroups - 10 > ndpts / 2 ) {
          return -1;
      }

      /* Early test in particular case for more general test below */
      /* "Test to see if the group widths and lengths are consistent with number of */
      /*  values, and length of section 7. */
      if( idrstmpl[12] < 0 || idrstmpl[14] < 0 || idrstmpl[14] > ndpts )
          return -1;
      if( nbitsglen == 0 &&
          ((ngroups > 1 && idrstmpl[12] != (ndpts - idrstmpl[14]) / (ngroups - 1)) ||
           idrstmpl[12] * (ngroups-1) + idrstmpl[14] != ndpts) )
      {
          return -1;
      }

      iofst=0;
      ifld=(g2int *)calloc(ndpts,sizeof(g2int));
      //printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
      gref=(g2int *)calloc(ngroups,sizeof(g2int));
      //printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
      gwidth=(g2int *)calloc(ngroups,sizeof(g2int));
      //printf("ALLOC gwidth: %d %x\n",(int)ngroups,gwidth);
      if( ifld == NULL || gref == NULL || gwidth == NULL )
      {
          free(ifld);
          free(gref);
          free(gwidth);
          return -1;
      }
//
//  Get missing values, if supplied
//
      if ( idrstmpl[6] == 1 ) {
         if (itype == 0)
            rdieee(idrstmpl+7,&rmiss1,1);
         else
            rmiss1=(g2float)idrstmpl[7];
      }
      if ( idrstmpl[6] == 2 ) {
         if (itype == 0) {
            rdieee(idrstmpl+7,&rmiss1,1);
            rdieee(idrstmpl+8,&rmiss2,1);
         }
         else {
            rmiss1=(g2float)idrstmpl[7];
            rmiss2=(g2float)idrstmpl[8];
         }
      }

      //printf("RMISSs: %f %f %f \n",rmiss1,rmiss2,ref);
//
//  Extract spatial differencing values, if using DRS Template 5.3
//
      if (idrsnum == 3) {
         if (nbitsd != 0) {
// one mistake here should be unsigned int
              int ok = gbit2(cpack, cpack_length, &ival1, iofst,nbitsd) != -1;
              iofst=iofst+nbitsd;
//              gbit(cpack,&isign,iofst,1);
//              iofst=iofst+1
//              gbit(cpack,&ival1,iofst,nbitsd-1);
//              iofst=iofst+nbitsd-1;
//              if (isign == 1) ival1=-ival1;
              if (idrstmpl[16] == 2) {
// one mistake here should be unsigned int
                 ok = ok && gbit2(cpack, cpack_length, &ival2,iofst,nbitsd) != -1;
                 iofst=iofst+nbitsd;
//                 gbit(cpack,&isign,iofst,1);
//                 iofst=iofst+1;
//                 gbit(cpack,&ival2,iofst,nbitsd-1);
//                 iofst=iofst+nbitsd-1;
//                 if (isign == 1) ival2=-ival2;
              }
              ok = ok && gbit2(cpack, cpack_length, &isign,iofst,1) != -1;
              iofst=iofst+1;
              ok = ok && gbit2(cpack, cpack_length, &minsd,iofst,nbitsd-1) != - 1;
              iofst=iofst+nbitsd-1;

              if (!ok)
              {
                  free(ifld);
                  free(gref);
                  free(gwidth);
                  return -1;
              }

              if (isign == 1 && minsd != INT_MIN) minsd=-minsd;
         }
         else {
              ival1=0;
              ival2=0;
              minsd=0;
         }
       //printf("SDu %ld %ld %ld %ld \n",ival1,ival2,minsd,nbitsd);
      }
//
//  Extract Each Group's reference value
//
      //printf("SAG1: %ld %ld %ld \n",nbitsgref,ngroups,iofst);
      if (nbitsgref != 0) {
         if( gbits(cpack,cpack_length,gref+0,iofst,nbitsgref,0,ngroups) != 0 )
         {
             free(ifld);
             free(gwidth);
             free(gref);
             return -1;
         }
         itemp=nbitsgref*ngroups;
         iofst=iofst+itemp;
         if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
      }
#if 0
      else {
         for (j=0;j<ngroups;j++)
              gref[j]=0;
      }
#endif

//
//  Extract Each Group's bit width
//
      //printf("SAG2: %ld %ld %ld %ld \n",nbitsgwidth,ngroups,iofst,idrstmpl[10]);
      if (nbitsgwidth != 0) {
         if( gbits(cpack,cpack_length,gwidth+0,iofst,nbitsgwidth,0,ngroups) != 0 )
         {
             free(ifld);
             free(gwidth);
             free(gref);
             return -1;
         }
         itemp=nbitsgwidth*ngroups;
         iofst=iofst+itemp;
         if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
      }
#if 0
      else {
         for (j=0;j<ngroups;j++)
                gwidth[j]=0;
      }
#endif

      for (j=0;j<ngroups;j++)
      {
          if( gwidth[j] > INT_MAX - idrstmpl[10] )
          {
             free(ifld);
             free(gwidth);
             free(gref);
             return -1;
          }
          gwidth[j]=gwidth[j]+idrstmpl[10];
      }

//
//  Extract Each Group's length (number of values in each group)
//
      glen=(g2int *)calloc(ngroups,sizeof(g2int));
      if( glen == NULL )
      {
        free(ifld);
        free(gwidth);
        free(gref);
        return -1;
      }
      //printf("ALLOC glen: %d %x\n",(int)ngroups,glen);
      //printf("SAG3: %ld %ld %ld %ld %ld \n",nbitsglen,ngroups,iofst,idrstmpl[13],idrstmpl[12]);
      if (nbitsglen != 0) {
         if( gbits(cpack,cpack_length,glen,iofst,nbitsglen,0,ngroups) != 0 )
         {
            free(ifld);
            free(gwidth);
            free(glen);
            free(gref);
             return -1;
         }
         itemp=nbitsglen*ngroups;
         iofst=iofst+itemp;
         if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
      }
#if 0
      else {
         for (j=0;j<ngroups;j++)
              glen[j]=0;
      }
#endif

      // Note: iterate only to ngroups - 1 since we override just after the
      // loop. Otherwise the security checks in the loop might trigger, like
      // on band 23 of
      // https://nomads.ncep.noaa.gov/pub/data/nccf/com/blend/prod/blend.20200605/16/grib2/blend.t16z.master.f001.co.grib2
      for (j=0;j<ngroups - 1;j++)
      {
           if( glen[j] < 0 ||
               (idrstmpl[13] != 0 && glen[j] > INT_MAX / idrstmpl[13]) ||
               glen[j] *  idrstmpl[13] > INT_MAX - idrstmpl[12] )
           {
                free(ifld);
                free(gwidth);
                free(glen);
                free(gref);
                return -1;
           }
           glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
      }
      glen[ngroups-1]=idrstmpl[14];
//
//  Test to see if the group widths and lengths are consistent with number of
//  values, and length of section 7.
//
      totBit = 0;
      totLen = 0;
      for (j=0;j<ngroups;j++) {
        g2int width_mult_len;
        if( gwidth[j] < 0 || glen[j] < 0 ||
            (gwidth[j] > 0 && glen[j] > INT_MAX / gwidth[j]) )
        {
            errorOccurred = 1;
            break;
        }
        width_mult_len = gwidth[j]*glen[j];
        if( totBit > INT_MAX - width_mult_len )
        {
            errorOccurred = 1;
            break;
        }
        totBit += width_mult_len;
        if( totLen > INT_MAX - glen[j] )
        {
            errorOccurred = 1;
            break;
        }
        totLen += glen[j];
      }
      if (errorOccurred || totLen != ndpts || totBit / 8. > lensec) {
        free(ifld);
        free(gwidth);
        free(glen);
        free(gref);
        return 1;
      }
//
//  For each group, unpack data values
//
      if ( idrstmpl[6] == 0 ) {        // no missing values
         n=0;
         for (j=0;j<ngroups;j++) {
           if (gwidth[j] != 0) {
             if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
             {
                 free(ifld);
                 free(gwidth);
                 free(glen);
                 free(gref);
                 return -1;
             }
             for (k=0;k<glen[j];k++) {
               ifld[n]=ifld[n]+gref[j];
               n=n+1;
             }
           }
           else {
             for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
             n=n+glen[j];
           }
           iofst=iofst+(gwidth[j]*glen[j]);
         }
      }
      else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
         // missing values included
         ifldmiss=(g2int *)malloc(ndpts*sizeof(g2int));
         //printf("ALLOC ifldmiss: %d %x\n",(int)ndpts,ifldmiss);
         //for (j=0;j<ndpts;j++) ifldmiss[j]=0;
         n=0;
         non=0;
         for (j=0;j<ngroups;j++) {
           //printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
           if (gwidth[j] != 0) {
             msng1=(g2int)int_power(2.0,gwidth[j])-1;
             msng2=msng1-1;
             if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
             {
                 free(ifld);
                 free(gwidth);
                 free(glen);
                 free(gref);
                 free(ifldmiss);
                 return -1;
             }
             iofst=iofst+(gwidth[j]*glen[j]);
             for (k=0;k<glen[j];k++) {
               if (ifld[n] == msng1) {
                  ifldmiss[n]=1;
                  //ifld[n]=0;
               }
               else if (idrstmpl[6]==2 && ifld[n]==msng2) {
                  ifldmiss[n]=2;
                  //ifld[n]=0;
               }
               else {
                  ifldmiss[n]=0;
                  ifld[non++]=ifld[n]+gref[j];
               }
               n++;
             }
           }
           else {
             msng1=(g2int)int_power(2.0,nbitsgref)-1;
             msng2=msng1-1;
             if (gref[j] == msng1) {
                for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
             }
             else if (idrstmpl[6]==2 && gref[j]==msng2) {
                for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
             }
             else {
                for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
                for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
                non += glen[j];
             }
             n=n+glen[j];
           }
         }
      }

      if ( gref != 0 ) free(gref);
      if ( gwidth != 0 ) free(gwidth);
      if ( glen != 0 ) free(glen);
//
//  If using spatial differences, add overall min value, and
//  sum up recursively
//
      //printf("SAGod: %ld %ld\n",idrsnum,idrstmpl[16]);
      if (idrsnum == 3) {         // spatial differencing
         if (idrstmpl[16] == 1) {      // first order
            ifld[0]=ival1;
            if ( idrstmpl[6] == 0 ) itemp=ndpts;        // no missing values
            else  itemp=non;
            for (n=1;n<itemp;n++) {
               if( (minsd > 0 && ifld[n] > INT_MAX - minsd) ||
                   (minsd < 0 && ifld[n] < INT_MIN - minsd) )
               {
                   free(ifldmiss);
                   free(ifld);
                   return -1;
               }
               ifld[n]=ifld[n]+minsd;
               if( (ifld[n-1] > 0 && ifld[n] > INT_MAX - ifld[n-1]) ||
                   (ifld[n-1] < 0 && ifld[n] < INT_MIN - ifld[n-1]) )
               {
                   free(ifldmiss);
                   free(ifld);
                   return -1;
               }
               ifld[n]=ifld[n]+ifld[n-1];
            }
         }
         else if (idrstmpl[16] == 2) {    // second order
            ifld[0]=ival1;
            ifld[1]=ival2;
            if ( idrstmpl[6] == 0 ) itemp=ndpts;        // no missing values
            else  itemp=non;
            for (n=2;n<itemp;n++) {
               /* Lazy way of detecting int overflows: operate on double */
               double tmp = (double)ifld[n]+(double)minsd+(2.0 * ifld[n-1])-ifld[n-2];
               if( tmp > INT_MAX || tmp < INT_MIN )
               {
                   free(ifldmiss);
                   free(ifld);
                   return -1;
               }
               ifld[n]=(int)tmp;
            }
         }
      }
//
//  Scale data back to original form
//
      //printf("SAGT: %f %f %f\n",ref,bscale,dscale);
      if ( idrstmpl[6] == 0 ) {        // no missing values
         for (n=0;n<ndpts;n++) {
            fld[n]=(((g2float)ifld[n]*bscale)+ref)*dscale;
         }
      }
      else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
         // missing values included
         non=0;
         for (n=0;n<ndpts;n++) {
            if ( ifldmiss[n] == 0 ) {
               fld[n]=(((g2float)ifld[non++]*bscale)+ref)*dscale;
               //printf(" SAG %d %f %d %f %f %f\n",n,fld[n],ifld[non-1],bscale,ref,dscale);
            }
            else if ( ifldmiss[n] == 1 )
               fld[n]=rmiss1;
            else if ( ifldmiss[n] == 2 )
               fld[n]=rmiss2;
         }
         if ( ifldmiss != 0 ) free(ifldmiss);
      }

      if ( ifld != 0 ) free(ifld);

      return(0);

}

Coverage Report

Created: 2026-09-26 08:22

Line	Count	Source
1		#include <float.h>
2		#include <limits.h>
3		#include <stdio.h>
4		#include <stdlib.h>
5		#include "grib2.h"
6
7		#ifndef DoubleToFloatClamp_defined
8		#define DoubleToFloatClamp_defined
9	818	static float DoubleToFloatClamp(double val) {
10	818	if (val >= FLT_MAX) return FLT_MAX;
11	818	if (val <= -FLT_MAX) return -FLT_MAX;
12	818	return (float)val;
13	818	}
14		#endif
15
16		int comunpack(unsigned char cpack,g2int cpack_length,g2int lensec,g2int idrsnum,g2int idrstmpl,g2int ndpts,g2float *fld)
17		////$$$ SUBPROGRAM DOCUMENTATION BLOCK
18		// . . . .
19		// SUBPROGRAM: comunpack
20		// PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-10-29
21		//
22		// ABSTRACT: This subroutine unpacks a data field that was packed using a
23		// complex packing algorithm as defined in the GRIB2 documentation,
24		// using info from the GRIB2 Data Representation Template 5.2 or 5.3.
25		// Supports GRIB2 complex packing templates with or without
26		// spatial differences (i.e. DRTs 5.2 and 5.3).
27		//
28		// PROGRAM HISTORY LOG:
29		// 2002-10-29 Gilbert
30		// 2004-12-16 Gilbert - Added test ( provided by Arthur Taylor/MDL )
31		// to verify that group widths and lengths are
32		// consistent with section length.
33		//
34		// USAGE: int comunpack(unsigned char *cpack,g2int lensec,g2int idrsnum,
35		// g2int idrstmpl, g2int ndpts,g2float fld)
36		// INPUT ARGUMENT LIST:
37		// cpack - pointer to the packed data field.
38		// lensec - length of section 7 (used for error checking).
39		// idrsnum - Data Representation Template number 5.N
40		// Must equal 2 or 3.
41		// idrstmpl - pointer to the array of values for Data Representation
42		// Template 5.2 or 5.3
43		// ndpts - The number of data values to unpack
44		//
45		// OUTPUT ARGUMENT LIST:
46		// fld - Contains the unpacked data values. fld must be allocated
47		// with at least ndpts*sizeof(g2float) bytes before
48		// calling this routine.
49		//
50		// REMARKS: None
51		//
52		// ATTRIBUTES:
53		// LANGUAGE: C
54		// MACHINE:
55		//
56		//$$$//
57	409	{
58
59	409	g2int nbitsd=0,isign;
60	409	g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
61	409	g2int ifld=NULL,ifldmiss=NULL;
62	409	g2int gref=NULL,gwidth=NULL,*glen=NULL;
63	409	g2int itype,ngroups,nbitsgref,nbitsgwidth,nbitsglen;
64	409	g2int msng1,msng2;
65	409	g2float ref,bscale,dscale,rmiss1,rmiss2;
66	409	g2int totBit, totLen;
67	409	g2int errorOccurred = 0;
68
69		//printf('IDRSTMPL: ',(idrstmpl(j),j=1,16)
70	409	rdieee(idrstmpl+0,&ref,1);
71		// printf("SAGTref: %f\n",ref);
72	409	bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
73	409	dscale = DoubleToFloatClamp(int_power(10.0,-idrstmpl[2]));
74	409	nbitsgref = idrstmpl[3];
75	409	itype = idrstmpl[4];
76	409	ngroups = idrstmpl[9];
77	409	nbitsgwidth = idrstmpl[11];
78	409	nbitsglen = idrstmpl[15];
79	409	if (idrsnum == 3)
80	77	nbitsd=idrstmpl[17]*8;
81
82		// Constant field
83
84	409	if (ngroups == 0) {
85	0	for (j=0;j<ndpts;j++) fld[j]=ref;
86	0	return(0);
87	0	}
88
89		/* To avoid excessive memory allocations. Not completely sure */
90		/* if this test is appropriate (the 10 and 2 are arbitrary), */
91		/* but it doesn't seem to make sense to have ngroups much larger than */
92		/* ndpts */
93	409	if( ngroups < 0 \|\| ngroups - 10 > ndpts / 2 ) {
94	27	return -1;
95	27	}
96
97		/* Early test in particular case for more general test below */
98		/* "Test to see if the group widths and lengths are consistent with number of */
99		/* values, and length of section 7. */
100	382	if( idrstmpl[12] < 0 \|\| idrstmpl[14] < 0 \|\| idrstmpl[14] > ndpts )
101	10	return -1;
102	372	if( nbitsglen == 0 &&
103	0	((ngroups > 1 && idrstmpl[12] != (ndpts - idrstmpl[14]) / (ngroups - 1)) \|\|
104	0	idrstmpl[12] * (ngroups-1) + idrstmpl[14] != ndpts) )
105	0	{
106	0	return -1;
107	0	}
108
109	372	iofst=0;
110	372	ifld=(g2int *)calloc(ndpts,sizeof(g2int));
111		//printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
112	372	gref=(g2int *)calloc(ngroups,sizeof(g2int));
113		//printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
114	372	gwidth=(g2int *)calloc(ngroups,sizeof(g2int));
115		//printf("ALLOC gwidth: %d %x\n",(int)ngroups,gwidth);
116	372	if( ifld == NULL \|\| gref == NULL \|\| gwidth == NULL )
117	0	{
118	0	free(ifld);
119	0	free(gref);
120	0	free(gwidth);
121	0	return -1;
122	0	}
123		//
124		// Get missing values, if supplied
125		//
126	372	if ( idrstmpl[6] == 1 ) {
127	341	if (itype == 0)
128	341	rdieee(idrstmpl+7,&rmiss1,1);
129	0	else
130	0	rmiss1=(g2float)idrstmpl[7];
131	341	}
132	372	if ( idrstmpl[6] == 2 ) {
133	0	if (itype == 0) {
134	0	rdieee(idrstmpl+7,&rmiss1,1);
135	0	rdieee(idrstmpl+8,&rmiss2,1);
136	0	}
137	0	else {
138	0	rmiss1=(g2float)idrstmpl[7];
139	0	rmiss2=(g2float)idrstmpl[8];
140	0	}
141	0	}
142
143		//printf("RMISSs: %f %f %f \n",rmiss1,rmiss2,ref);
144		//
145		// Extract spatial differencing values, if using DRS Template 5.3
146		//
147	372	if (idrsnum == 3) {
148	76	if (nbitsd != 0) {
149		// one mistake here should be unsigned int
150	76	int ok = gbit2(cpack, cpack_length, &ival1, iofst,nbitsd) != -1;
151	76	iofst=iofst+nbitsd;
152		// gbit(cpack,&isign,iofst,1);
153		// iofst=iofst+1
154		// gbit(cpack,&ival1,iofst,nbitsd-1);
155		// iofst=iofst+nbitsd-1;
156		// if (isign == 1) ival1=-ival1;
157	76	if (idrstmpl[16] == 2) {
158		// one mistake here should be unsigned int
159	50	ok = ok && gbit2(cpack, cpack_length, &ival2,iofst,nbitsd) != -1;
160	50	iofst=iofst+nbitsd;
161		// gbit(cpack,&isign,iofst,1);
162		// iofst=iofst+1;
163		// gbit(cpack,&ival2,iofst,nbitsd-1);
164		// iofst=iofst+nbitsd-1;
165		// if (isign == 1) ival2=-ival2;
166	50	}
167	76	ok = ok && gbit2(cpack, cpack_length, &isign,iofst,1) != -1;
168	76	iofst=iofst+1;
169	76	ok = ok && gbit2(cpack, cpack_length, &minsd,iofst,nbitsd-1) != - 1;
170	76	iofst=iofst+nbitsd-1;
171
172	76	if (!ok)
173	0	{
174	0	free(ifld);
175	0	free(gref);
176	0	free(gwidth);
177	0	return -1;
178	0	}
179
180	76	if (isign == 1 && minsd != INT_MIN) minsd=-minsd;
181	76	}
182	0	else {
183	0	ival1=0;
184	0	ival2=0;
185	0	minsd=0;
186	0	}
187		//printf("SDu %ld %ld %ld %ld \n",ival1,ival2,minsd,nbitsd);
188	76	}
189		//
190		// Extract Each Group's reference value
191		//
192		//printf("SAG1: %ld %ld %ld \n",nbitsgref,ngroups,iofst);
193	372	if (nbitsgref != 0) {
194	372	if( gbits(cpack,cpack_length,gref+0,iofst,nbitsgref,0,ngroups) != 0 )
195	10	{
196	10	free(ifld);
197	10	free(gwidth);
198	10	free(gref);
199	10	return -1;
200	10	}
201	362	itemp=nbitsgref*ngroups;
202	362	iofst=iofst+itemp;
203	362	if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
204	362	}
205		#if 0
206		else {
207		for (j=0;j<ngroups;j++)
208		gref[j]=0;
209		}
210		#endif
211
212		//
213		// Extract Each Group's bit width
214		//
215		//printf("SAG2: %ld %ld %ld %ld \n",nbitsgwidth,ngroups,iofst,idrstmpl[10]);
216	362	if (nbitsgwidth != 0) {
217	362	if( gbits(cpack,cpack_length,gwidth+0,iofst,nbitsgwidth,0,ngroups) != 0 )
218	0	{
219	0	free(ifld);
220	0	free(gwidth);
221	0	free(gref);
222	0	return -1;
223	0	}
224	362	itemp=nbitsgwidth*ngroups;
225	362	iofst=iofst+itemp;
226	362	if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
227	362	}
228		#if 0
229		else {
230		for (j=0;j<ngroups;j++)
231		gwidth[j]=0;
232		}
233		#endif
234
235	778k	for (j=0;j<ngroups;j++)
236	778k	{
237	778k	if( gwidth[j] > INT_MAX - idrstmpl[10] )
238	0	{
239	0	free(ifld);
240	0	free(gwidth);
241	0	free(gref);
242	0	return -1;
243	0	}
244	778k	gwidth[j]=gwidth[j]+idrstmpl[10];
245	778k	}
246
247		//
248		// Extract Each Group's length (number of values in each group)
249		//
250	362	glen=(g2int *)calloc(ngroups,sizeof(g2int));
251	362	if( glen == NULL )
252	0	{
253	0	free(ifld);
254	0	free(gwidth);
255	0	free(gref);
256	0	return -1;
257	0	}
258		//printf("ALLOC glen: %d %x\n",(int)ngroups,glen);
259		//printf("SAG3: %ld %ld %ld %ld %ld \n",nbitsglen,ngroups,iofst,idrstmpl[13],idrstmpl[12]);
260	362	if (nbitsglen != 0) {
261	362	if( gbits(cpack,cpack_length,glen,iofst,nbitsglen,0,ngroups) != 0 )
262	16	{
263	16	free(ifld);
264	16	free(gwidth);
265	16	free(glen);
266	16	free(gref);
267	16	return -1;
268	16	}
269	346	itemp=nbitsglen*ngroups;
270	346	iofst=iofst+itemp;
271	346	if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
272	346	}
273		#if 0
274		else {
275		for (j=0;j<ngroups;j++)
276		glen[j]=0;
277		}
278		#endif
279
280		// Note: iterate only to ngroups - 1 since we override just after the
281		// loop. Otherwise the security checks in the loop might trigger, like
282		// on band 23 of
283		// https://nomads.ncep.noaa.gov/pub/data/nccf/com/blend/prod/blend.20200605/16/grib2/blend.t16z.master.f001.co.grib2
284	777k	for (j=0;j<ngroups - 1;j++)
285	776k	{
286	776k	if( glen[j] < 0 \|\|
287	776k	(idrstmpl[13] != 0 && glen[j] > INT_MAX / idrstmpl[13]) \|\|
288	776k	glen[j] * idrstmpl[13] > INT_MAX - idrstmpl[12] )
289	0	{
290	0	free(ifld);
291	0	free(gwidth);
292	0	free(glen);
293	0	free(gref);
294	0	return -1;
295	0	}
296	776k	glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
297	776k	}
298	346	glen[ngroups-1]=idrstmpl[14];
299		//
300		// Test to see if the group widths and lengths are consistent with number of
301		// values, and length of section 7.
302		//
303	346	totBit = 0;
304	346	totLen = 0;
305	777k	for (j=0;j<ngroups;j++) {
306	777k	g2int width_mult_len;
307	777k	if( gwidth[j] < 0 \|\| glen[j] < 0 \|\|
308	777k	(gwidth[j] > 0 && glen[j] > INT_MAX / gwidth[j]) )
309	0	{
310	0	errorOccurred = 1;
311	0	break;
312	0	}
313	777k	width_mult_len = gwidth[j]*glen[j];
314	777k	if( totBit > INT_MAX - width_mult_len )
315	0	{
316	0	errorOccurred = 1;
317	0	break;
318	0	}
319	777k	totBit += width_mult_len;
320	777k	if( totLen > INT_MAX - glen[j] )
321	0	{
322	0	errorOccurred = 1;
323	0	break;
324	0	}
325	777k	totLen += glen[j];
326	777k	}
327	346	if (errorOccurred \|\| totLen != ndpts \|\| totBit / 8. > lensec) {
328	105	free(ifld);
329	105	free(gwidth);
330	105	free(glen);
331	105	free(gref);
332	105	return 1;
333	105	}
334		//
335		// For each group, unpack data values
336		//
337	241	if ( idrstmpl[6] == 0 ) { // no missing values
338	20	n=0;
339	254k	for (j=0;j<ngroups;j++) {
340	254k	if (gwidth[j] != 0) {
341	239k	if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
342	0	{
343	0	free(ifld);
344	0	free(gwidth);
345	0	free(glen);
346	0	free(gref);
347	0	return -1;
348	0	}
349	3.73M	for (k=0;k<glen[j];k++) {
350	3.49M	ifld[n]=ifld[n]+gref[j];
351	3.49M	n=n+1;
352	3.49M	}
353	239k	}
354	15.4k	else {
355	33.4k	for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
356	15.4k	n=n+glen[j];
357	15.4k	}
358	254k	iofst=iofst+(gwidth[j]*glen[j]);
359	254k	}
360	20	}
361	221	else if ( idrstmpl[6]==1 \|\| idrstmpl[6]==2 ) {
362		// missing values included
363	220	ifldmiss=(g2int )malloc(ndptssizeof(g2int));
364		//printf("ALLOC ifldmiss: %d %x\n",(int)ndpts,ifldmiss);
365		//for (j=0;j<ndpts;j++) ifldmiss[j]=0;
366	220	n=0;
367	220	non=0;
368	512k	for (j=0;j<ngroups;j++) {
369		//printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
370	512k	if (gwidth[j] != 0) {
371	446k	msng1=(g2int)int_power(2.0,gwidth[j])-1;
372	446k	msng2=msng1-1;
373	446k	if( gbits(cpack,cpack_length,ifld+n,iofst,gwidth[j],0,glen[j]) != 0 )
374	8	{
375	8	free(ifld);
376	8	free(gwidth);
377	8	free(glen);
378	8	free(gref);
379	8	free(ifldmiss);
380	8	return -1;
381	8	}
382	446k	iofst=iofst+(gwidth[j]*glen[j]);
383	6.44M	for (k=0;k<glen[j];k++) {
384	5.99M	if (ifld[n] == msng1) {
385	576k	ifldmiss[n]=1;
386		//ifld[n]=0;
387	576k	}
388	5.42M	else if (idrstmpl[6]==2 && ifld[n]==msng2) {
389	0	ifldmiss[n]=2;
390		//ifld[n]=0;
391	0	}
392	5.42M	else {
393	5.42M	ifldmiss[n]=0;
394	5.42M	ifld[non++]=ifld[n]+gref[j];
395	5.42M	}
396	5.99M	n++;
397	5.99M	}
398	446k	}
399	65.2k	else {
400	65.2k	msng1=(g2int)int_power(2.0,nbitsgref)-1;
401	65.2k	msng2=msng1-1;
402	65.2k	if (gref[j] == msng1) {
403	7.12M	for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
404	64.8k	}
405	384	else if (idrstmpl[6]==2 && gref[j]==msng2) {
406	0	for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
407	0	}
408	384	else {
409	5.33k	for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
410	5.33k	for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
411	384	non += glen[j];
412	384	}
413	65.2k	n=n+glen[j];
414	65.2k	}
415	512k	}
416	220	}
417
418	233	if ( gref != 0 ) free(gref);
419	233	if ( gwidth != 0 ) free(gwidth);
420	233	if ( glen != 0 ) free(glen);
421		//
422		// If using spatial differences, add overall min value, and
423		// sum up recursively
424		//
425		//printf("SAGod: %ld %ld\n",idrsnum,idrstmpl[16]);
426	233	if (idrsnum == 3) { // spatial differencing
427	66	if (idrstmpl[16] == 1) { // first order
428	24	ifld[0]=ival1;
429	24	if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
430	24	else itemp=non;
431	442k	for (n=1;n<itemp;n++) {
432	442k	if( (minsd > 0 && ifld[n] > INT_MAX - minsd) \|\|
433	442k	(minsd < 0 && ifld[n] < INT_MIN - minsd) )
434	0	{
435	0	free(ifldmiss);
436	0	free(ifld);
437	0	return -1;
438	0	}
439	442k	ifld[n]=ifld[n]+minsd;
440	442k	if( (ifld[n-1] > 0 && ifld[n] > INT_MAX - ifld[n-1]) \|\|
441	442k	(ifld[n-1] < 0 && ifld[n] < INT_MIN - ifld[n-1]) )
442	0	{
443	0	free(ifldmiss);
444	0	free(ifld);
445	0	return -1;
446	0	}
447	442k	ifld[n]=ifld[n]+ifld[n-1];
448	442k	}
449	24	}
450	42	else if (idrstmpl[16] == 2) { // second order
451	42	ifld[0]=ival1;
452	42	ifld[1]=ival2;
453	42	if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
454	28	else itemp=non;
455	3.12M	for (n=2;n<itemp;n++) {
456		/* Lazy way of detecting int overflows: operate on double */
457	3.12M	double tmp = (double)ifld[n]+(double)minsd+(2.0 * ifld[n-1])-ifld[n-2];
458	3.12M	if( tmp > INT_MAX \|\| tmp < INT_MIN )
459	8	{
460	8	free(ifldmiss);
461	8	free(ifld);
462	8	return -1;
463	8	}
464	3.12M	ifld[n]=(int)tmp;
465	3.12M	}
466	42	}
467	66	}
468		//
469		// Scale data back to original form
470		//
471		//printf("SAGT: %f %f %f\n",ref,bscale,dscale);
472	225	if ( idrstmpl[6] == 0 ) { // no missing values
473	1.50M	for (n=0;n<ndpts;n++) {
474	1.50M	fld[n]=(((g2float)ifld[n]bscale)+ref)dscale;
475	1.50M	}
476	12	}
477	213	else if ( idrstmpl[6]==1 \|\| idrstmpl[6]==2 ) {
478		// missing values included
479	212	non=0;
480	5.38M	for (n=0;n<ndpts;n++) {
481	5.38M	if ( ifldmiss[n] == 0 ) {
482	2.45M	fld[n]=(((g2float)ifld[non++]bscale)+ref)dscale;
483		//printf(" SAG %d %f %d %f %f %f\n",n,fld[n],ifld[non-1],bscale,ref,dscale);
484	2.45M	}
485	2.92M	else if ( ifldmiss[n] == 1 )
486	2.92M	fld[n]=rmiss1;
487	0	else if ( ifldmiss[n] == 2 )
488	0	fld[n]=rmiss2;
489	5.38M	}
490	212	if ( ifldmiss != 0 ) free(ifldmiss);
491	212	}
492
493	225	if ( ifld != 0 ) free(ifld);
494
495	225	return(0);
496
497	233	}