grib/src/main/java/ucar/nc2/grib/grib2/Grib2DataReader.java
private float[] getData3(RandomAccessFile raf, Grib2Drs.Type3 gdrs) throws IOException {
    int mvm = gdrs.missingValueManagement;
    float mv = getMissingValue(gdrs);

    float DD = (float) java.lang.Math.pow((double) 10, (double) gdrs.decimalScaleFactor);
    float R = gdrs.referenceValue;
    float EE = (float) java.lang.Math.pow(2.0, (double) gdrs.binaryScaleFactor);
    float ref_val = R / DD;

    int NG = gdrs.numberOfGroups;
    if (NG == 0) {
      return nGroups0(bitmapIndicator, ref_val, mv);
    }

    BitReader reader = new BitReader(raf, startPos + 5);

    int ival1;
    int ival2 = 0;
    int minsd;

    // [6-ww]   1st values of undifferenced scaled values and minimums
    int os = gdrs.orderSpatial;
    int nbitsd = gdrs.descriptorSpatial;
    int sign;
    // ds is number of bytes, convert to bits -1 for sign bit
    nbitsd = nbitsd * 8;
    if (nbitsd > 0) {         // first order spatial differencing g1 and gMin
      sign = (int) reader.bits2UInt(1);
      ival1 = (int) reader.bits2UInt(nbitsd - 1);
      if (sign == 1) {
        ival1 = -ival1;
      }
      if (os == 2) {  //second order spatial differencing h1, h2, hMin
        sign = (int) reader.bits2UInt(1);
        ival2 = (int) reader.bits2UInt(nbitsd - 1);
        if (sign == 1) {
          ival2 = -ival2;
        }
      }
      sign = (int) reader.bits2UInt(1);
      minsd = (int) reader.bits2UInt(nbitsd - 1);
      if (sign == 1) {
        minsd = -minsd;
      }

    } else {
      float[] data = new float[totalNPoints];
      for (int i = 0; i < totalNPoints; i++) {
        data[i] = mv;
      }
      return data;
    }

    // [ww +1]-xx  Get reference values for groups (X1's)
    // X1 == gref
    int[] X1 = new int[NG]; // initialized to zero
    int nb = gdrs.numberOfBits;
    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        X1[i] = (int) reader.bits2UInt(nb);
      }
    }

    // [xx +1 ]-yy Get number of bits used to encode each group
    // NB == gwidth
    int[] NB = new int[NG]; // initialized to zero
    nb = gdrs.bitsGroupWidths;
    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        NB[i] = (int) reader.bits2UInt(nb);
      }
    }

    int referenceGroupWidths = gdrs.referenceGroupWidths;
    for (int i = 0; i < NG; i++) {
      NB[i] += referenceGroupWidths;
    }

    // [yy +1 ]-zz Get the scaled group lengths using formula
    //     Ln = ref + Kn * len_inc, where n = 1-NG,
    //          ref = referenceGroupLength, and  len_inc = lengthIncrement

    int[] L = new int[NG];   // initialized to zero
    int referenceGroupLength = gdrs.referenceGroupLength;
    nb = gdrs.bitsScaledGroupLength;
    int len_inc = gdrs.lengthIncrement;

    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        L[i] = (int) reader.bits2UInt(nb);
      }
    }

    int totalL = 0;
    for (int i = 0; i < NG; i++) {
      L[i] = L[i] * len_inc + referenceGroupLength;
      totalL += L[i];
    }
    totalL -= L[NG - 1];
    totalL += gdrs.lengthLastGroup;

    //enter Length of Last Group
    L[NG - 1] = gdrs.lengthLastGroup;

    // test
    if (mvm != 0) {
      if (totalL != totalNPoints) {
        logger.warn("NPoints != gds.nPts: " + totalL + "!=" + totalNPoints);
        float[] data = new float[totalNPoints];
        for (int i = 0; i < totalNPoints; i++) {
          data[i] = mv;
        }
        return data;
      }
    } else {
      if (totalL != dataNPoints) {
        logger.warn("NPoints != drs.nPts: " + totalL + "!=" + totalNPoints);
        float[] data = new float[totalNPoints];
        for (int i = 0; i < totalNPoints; i++) {
          data[i] = mv;
        }
        return data;
      }
    }

    float[] data = new float[totalNPoints];

    // [zz +1 ]-nn get X2 values and calculate the results Y using formula
//      formula used to create values,  Y * 10**D = R + (X1 + X2) * 2**E

//               Y = (R + (X1 + X2) * (2 ** E) ) / (10 ** D)]
//               WHERE:
//                     Y = THE VALUE WE ARE UNPACKING
//                     R = THE REFERENCE VALUE (FIRST ORDER MINIMA)
//                    X1 = THE PACKED VALUE
//                    X2 = THE SECOND ORDER MINIMA
//                     E = THE BINARY SCALE FACTOR
//                     D = THE DECIMAL SCALE FACTOR
    int count = 0;
    reader.incrByte();
    int dataSize = 0;
    boolean[] dataBitMap = null;
    if (mvm == 0) {
      for (int i = 0; i < NG; i++) {
        if (NB[i] != 0) {
          for (int j = 0; j < L[i]; j++) {
            data[count++] = (int) reader.bits2UInt(NB[i]) + X1[i];
          }
        } else {
          for (int j = 0; j < L[i]; j++) {
            data[count++] = X1[i];
          }
        }
      }  // end for i

    } else if (mvm == 1 || mvm == 2) {
      // don't add missing values into data but keep track of them in dataBitMap
      dataBitMap = new boolean[totalNPoints];
      dataSize = 0;
      for (int i = 0; i < NG; i++) {
        if (NB[i] != 0) {
          int msng1 = bitsmv1[NB[i]];
          int msng2 = msng1 - 1;
          for (int j = 0; j < L[i]; j++) {
            data[count] = (int) reader.bits2UInt(NB[i]);
            if (data[count] == msng1 || mvm == 2 && data[count] == msng2) {
              dataBitMap[count] = false;
            } else {
              dataBitMap[count] = true;
              data[dataSize++] = data[count] + X1[i];
            }
            count++;
          }
        } else {  // (NB[i] == 0
          int msng1 = bitsmv1[gdrs.numberOfBits];
          int msng2 = msng1 - 1;
          if (X1[i] == msng1) {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count++] = false;
            }
            //data[count++] = X1[i];
          } else if (mvm == 2 && X1[i] == msng2) {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count++] = false;
            }
          } else {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count] = true;
              data[dataSize++] = X1[i];
              count++;
            }
          }
        }
      }  // end for i
    }

    // first order spatial differencing
    if (os == 1) {   // g1 and gMin
      // encoded by G(n) = F(n) - F(n -1 )
      // decoded by F(n) = G(n) + F(n -1 )
      // data[] at this point contains G0, G1, G2, ....
      data[0] = ival1;
      int itemp;
      if (mvm == 0) {           // no missing values
        itemp = totalNPoints;
      } else {
        itemp = dataSize;
      }
      for (int i = 1; i < itemp; i++) {
        data[i] += minsd;
        data[i] = data[i] + data[i - 1];
      }
    } else if (os == 2) { // 2nd order
      data[0] = ival1;
      data[1] = ival2;
      int itemp;
      if (mvm == 0) {           // no missing values
        itemp = totalNPoints;
      } else {
        itemp = dataSize;
      }
      for (int i = 2; i < itemp; i++) {
        data[i] += minsd;
        data[i] = data[i] + (2 * data[i - 1]) - data[i - 2];
      }
    }

    // formula used to create values,  Y * 10**D = R + (X1 + X2) * 2**E

    //               Y = (R + (X1 + X2) * (2 ** E) ) / (10 ** D)]
//               WHERE:
//                     Y = THE VALUE WE ARE UNPACKING
//                     R = THE REFERENCE VALUE (FIRST ORDER MINIMA)
//                    X1 = THE PACKED VALUE
//                    X2 = THE SECOND ORDER MINIMA
//                     E = THE BINARY SCALE FACTOR
//                     D = THE DECIMAL SCALE FACTOR

    if (mvm == 0) {  // no missing values
      for (int i = 0; i < data.length; i++) {
        data[i] = (R + (data[i] * EE)) / DD;
      }
    } else if (mvm == 1 || mvm == 2) {         // missing value == 1  || missing value == 2
      int count2 = 0;
      float[] tmp = new float[totalNPoints];
      for (int i = 0; i < data.length; i++) {
        if (dataBitMap[i]) {
          tmp[i] = (R + (data[count2++] * EE)) / DD;
        } else { // mvm = 1 or 2
          tmp[i] = mv;
        }
      }
      data = tmp;
    }

    // bit map is used
    if (bitmap != null) {
      int idx = 0;
      float[] tmp = new float[totalNPoints];
      for (int i = 0; i < totalNPoints; i++) {
        if (GribNumbers.testBitIsSet(bitmap[i / 8], i % 8)) {
          tmp[i] = data[idx++];
        } else {
          tmp[i] = mv;
        }
      }
      data = tmp;
    }

    return data;
  }
private float[] getData3(RandomAccessFile raf, Grib2Drs.Type3 gdrs) throws IOException {
    int mvm = gdrs.missingValueManagement;
    float mv = getMissingValue(gdrs);

    float DD = (float) java.lang.Math.pow((double) 10, (double) gdrs.decimalScaleFactor);
    float R = gdrs.referenceValue;
    float EE = (float) java.lang.Math.pow(2.0, (double) gdrs.binaryScaleFactor);
    float ref_val = R / DD;

    int NG = gdrs.numberOfGroups;
    if (NG == 0) {
      return nGroups0(bitmapIndicator, ref_val, mv);
    }

    BitReader reader = new BitReader(raf, startPos + 5);

    int ival1;
    int ival2 = 0;
    int minsd;

    // [6-ww]   1st values of undifferenced scaled values and minimums
    int os = gdrs.orderSpatial;
    int nbitsd = gdrs.descriptorSpatial;
    int sign;
    // ds is number of bytes, convert to bits -1 for sign bit
    nbitsd = nbitsd * 8;
    if (nbitsd > 0) {         // first order spatial differencing g1 and gMin
      sign = (int) reader.bits2UInt(1);
      ival1 = (int) reader.bits2UInt(nbitsd - 1);
      if (sign == 1) {
        ival1 = -ival1;
      }
      if (os == 2) {  //second order spatial differencing h1, h2, hMin
        sign = (int) reader.bits2UInt(1);
        ival2 = (int) reader.bits2UInt(nbitsd - 1);
        if (sign == 1) {
          ival2 = -ival2;
        }
      }
      sign = (int) reader.bits2UInt(1);
      minsd = (int) reader.bits2UInt(nbitsd - 1);
      if (sign == 1) {
        minsd = -minsd;
      }

    } else {
      float[] data = new float[totalNPoints];
      for (int i = 0; i < totalNPoints; i++) {
        data[i] = mv;
      }
      return data;
    }

    // [ww +1]-xx  Get reference values for groups (X1's)
    // X1 == gref
    int[] X1 = new int[NG]; // initialized to zero
    int nb = gdrs.numberOfBits;
    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        X1[i] = (int) reader.bits2UInt(nb);
      }
    }

    // [xx +1 ]-yy Get number of bits used to encode each group
    // NB == gwidth
    int[] NB = new int[NG]; // initialized to zero
    nb = gdrs.bitsGroupWidths;
    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        NB[i] = (int) reader.bits2UInt(nb);
      }
    }

    int referenceGroupWidths = gdrs.referenceGroupWidths;
    for (int i = 0; i < NG; i++) {
      NB[i] += referenceGroupWidths;
    }

    // [yy +1 ]-zz Get the scaled group lengths using formula
    //     Ln = ref + Kn * len_inc, where n = 1-NG,
    //          ref = referenceGroupLength, and  len_inc = lengthIncrement

    int[] L = new int[NG];   // initialized to zero
    int referenceGroupLength = gdrs.referenceGroupLength;
    nb = gdrs.bitsScaledGroupLength;
    int len_inc = gdrs.lengthIncrement;

    if (nb != 0) {
      reader.incrByte();
      for (int i = 0; i < NG; i++) {
        L[i] = (int) reader.bits2UInt(nb);
      }
    }

    int totalL = 0;
    for (int i = 0; i < NG; i++) {
      L[i] = L[i] * len_inc + referenceGroupLength;
      totalL += L[i];
    }
    totalL -= L[NG - 1];
    totalL += gdrs.lengthLastGroup;

    //enter Length of Last Group
    L[NG - 1] = gdrs.lengthLastGroup;

    // test
    if (mvm != 0) {
      if (totalL != totalNPoints) {
        logger.warn("NPoints != gds.nPts: " + totalL + "!=" + totalNPoints);
        float[] data = new float[totalNPoints];
        for (int i = 0; i < totalNPoints; i++) {
          data[i] = mv;
        }
        return data;
      }
    } else {
      if (totalL != dataNPoints) {
        logger.warn("NPoints != drs.nPts: " + totalL + "!=" + totalNPoints);
        float[] data = new float[totalNPoints];
        for (int i = 0; i < totalNPoints; i++) {
          data[i] = mv;
        }
        return data;
      }
    }

    float[] data = new float[totalNPoints];

    // [zz +1 ]-nn get X2 values and calculate the results Y using formula
//      formula used to create values,  Y * 10**D = R + (X1 + X2) * 2**E

//               Y = (R + (X1 + X2) * (2 ** E) ) / (10 ** D)]
//               WHERE:
//                     Y = THE VALUE WE ARE UNPACKING
//                     R = THE REFERENCE VALUE (FIRST ORDER MINIMA)
//                    X1 = THE PACKED VALUE
//                    X2 = THE SECOND ORDER MINIMA
//                     E = THE BINARY SCALE FACTOR
//                     D = THE DECIMAL SCALE FACTOR
    int count = 0;
    reader.incrByte();
    int dataSize = 0;
    boolean[] dataBitMap = null;
    if (mvm == 0) {
      for (int i = 0; i < NG; i++) {
        if (NB[i] != 0) {
          for (int j = 0; j < L[i]; j++) {
            data[count++] = (int) reader.bits2UInt(NB[i]) + X1[i];
          }
        } else {
          for (int j = 0; j < L[i]; j++) {
            data[count++] = X1[i];
          }
        }
      }  // end for i

    } else if (mvm == 1 || mvm == 2) {
      // don't add missing values into data but keep track of them in dataBitMap
      dataBitMap = new boolean[totalNPoints];
      dataSize = 0;
      for (int i = 0; i < NG; i++) {
        if (NB[i] != 0) {
          int msng1 = bitsmv1[NB[i]];
          int msng2 = msng1 - 1;
          for (int j = 0; j < L[i]; j++) {
            data[count] = (int) reader.bits2UInt(NB[i]);
            if (data[count] == msng1 || mvm == 2 && data[count] == msng2) {
              dataBitMap[count] = false;
            } else {
              dataBitMap[count] = true;
              data[dataSize++] = data[count] + X1[i];
            }
            count++;
          }
        } else {  // (NB[i] == 0
          int msng1 = bitsmv1[gdrs.numberOfBits];
          int msng2 = msng1 - 1;
          if (X1[i] == msng1) {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count++] = false;
            }
            //data[count++] = X1[i];
          } else if (mvm == 2 && X1[i] == msng2) {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count++] = false;
            }
          } else {
            for (int j = 0; j < L[i]; j++) {
              dataBitMap[count] = true;
              data[dataSize++] = X1[i];
              count++;
            }
          }
        }
      }  // end for i
    }

    // first order spatial differencing
    if (os == 1) {   // g1 and gMin
      // encoded by G(n) = F(n) - F(n -1 )
      // decoded by F(n) = G(n) + F(n -1 )
      // data[] at this point contains G0, G1, G2, ....
      data[0] = ival1;
      int itemp;
      if (mvm == 0) {           // no missing values
        itemp = totalNPoints;
      } else {
        itemp = dataSize;
      }
      for (int i = 1; i < itemp; i++) {
        data[i] += minsd;
        data[i] = data[i] + data[i - 1];
      }
    } else if (os == 2) { // 2nd order
      data[0] = ival1;
      data[1] = ival2;
      int itemp;
      if (mvm == 0) {           // no missing values
        itemp = totalNPoints;
      } else {
        itemp = dataSize;
      }
      for (int i = 2; i < itemp; i++) {
        data[i] += minsd;
        data[i] = data[i] + (2 * data[i - 1]) - data[i - 2];
      }
    }

    // formula used to create values,  Y * 10**D = R + (X1 + X2) * 2**E

    //               Y = (R + (X1 + X2) * (2 ** E) ) / (10 ** D)]
//               WHERE:
//                     Y = THE VALUE WE ARE UNPACKING
//                     R = THE REFERENCE VALUE (FIRST ORDER MINIMA)
//                    X1 = THE PACKED VALUE
//                    X2 = THE SECOND ORDER MINIMA
//                     E = THE BINARY SCALE FACTOR
//                     D = THE DECIMAL SCALE FACTOR

    if (mvm == 0) {  // no missing values
      for (int i = 0; i < data.length; i++) {
        data[i] = (R + (data[i] * EE)) / DD;
      }
    } else if (mvm == 1 || mvm == 2) {         // missing value == 1  || missing value == 2
      int count2 = 0;
      float[] tmp = new float[totalNPoints];
      for (int i = 0; i < data.length; i++) {
        if (dataBitMap[i]) {
          tmp[i] = (R + (data[count2++] * EE)) / DD;
        } else { // mvm = 1 or 2
          tmp[i] = mv;
        }
      }
      data = tmp;
    }

    // bit map is used
    if (bitmap != null) {
      int idx = 0;
      float[] tmp = new float[totalNPoints];
      for (int i = 0; i < totalNPoints; i++) {
        if (GribNumbers.testBitIsSet(bitmap[i / 8], i % 8)) {
          tmp[i] = data[idx++];
        } else {
          tmp[i] = mv;
        }
      }
      data = tmp;
    }

    return data;
  }
/*
Data template 7.3 – Grid point data – complex packing and spatial differencing
Note: For most templates, details of the packing process are described in Regulation 92.9.4.
Octet No. Contents
6–ww      First value(s) of original (undifferenced) scaled data values, followed by the overall minimum
of the differences. The number of values stored is 1 greater than the order of differentiation,
and the field width is described at octet 49 of data representation template 5.3 (see Note 1)
[ww+1]–xx NG group reference values (X1 in the decoding formula), each of which is encoded using the number
of bits specified in octet 20 of data representation template 5.0. Bits set to zero shall be
appended where necessary to ensure this sequence of numbers ends on an octet boundary
[xx+1]–nn Same as for data representation template 7.2

Notes:
(1) Referring to the notation in Note 1 of data representation template 5.3, at order 1, the values stored in octets 6–ww are
g1 and gmin. At order 2, the values stored are h1, h2, and hmin.
(2) Extra descriptors related to spatial differencing are added before the splitting descriptors, to reflect the separation
between the two approaches. It enables to share software parts between cases with and without spatial differencing.
(3) The position of overall minimum after initial data values is a choice that enables less software management.
(4) Overall minimum will be negative in most cases. First bit should indicate the sign: 0 if positive, 1 if negative.