Merge pull request #518 from NOAA-EMC/wgrib

wgrib src #445

Author: AlysonStahl-NOAA

wgrib/BDS_NValues.c

@@ -0,0 +1,28 @@
+#include <stdio.h>
+#include "bds.h"
+#include "grib.h"
+
+/*
+ * support for complex packing
+ *  determine the number of data points in the BDS
+ *  does not handle matrix values
+ */
+
+extern int ec_large_grib,  len_ec_bds;
+
+
+int BDS_NValues(unsigned char *bds) {
+
+    /* returns number of grid points as determined from the BDS */
+
+    int i = 0;
+
+    if (BDS_SimplePacking(bds) && BDS_Grid(bds)) {
+	i = ((BDS_LEN(bds) - BDS_DataStart(bds))*8 -
+		BDS_UnusedBits(bds)) / (BDS_NumBits(bds));
+    }
+    else if (BDS_ComplexPacking(bds) && BDS_Grid(bds)) {
+	i = BDS_P2(bds);
+    }
+    return i;
+}

wgrib/BDSunpk.c

@@ -0,0 +1,149 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stddef.h>
+#include <limits.h>
+#include "grib.h"
+#include "pds4.h"
+#include "bms.h"
+#include "bds.h"
+
+/* 1996				wesley ebisuzaki
+ *
+ * Unpack BDS section
+ *
+ * input: *bits, pointer to packed integer data
+ *        *bitmap, pointer to bitmap (undefined data), NULL if none
+ *        n_bits, number of bits per packed integer
+ *        n, number of data points (includes undefined data)
+ *        ref, scale: flt[] = ref + scale*packed_int
+ * output: *flt, pointer to output array
+ *        undefined values filled with UNDEFINED
+ *
+ * note: code assumes an integer > 32 bits
+ *
+ * 7/98 v1.2.1 fix bug for bitmaps and nbit >= 25 found by Larry Brasfield
+ * 2/01 v1.2.2 changed jj from long int to double
+ * 3/02 v1.2.3 added unpacking extensions for spectral data 
+ *             Luis Kornblueh, MPIfM 
+ * 7/06 v.1.2.4 fixed some bug complex packed data was not set to undefined
+ */
+
+static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
+static unsigned int map_masks[8] = {128, 64, 32, 16, 8, 4, 2, 1};
+static double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
+
+void BDS_unpack(float *flt, unsigned char *bds, unsigned char *bitmap,
+	int n_bits, int n, double ref, double scale) {
+
+    unsigned char *bits;
+
+    int i, mask_idx, t_bits, c_bits, j_bits;
+    unsigned int j, map_mask, tbits, jmask, bbits;
+    double jj;
+
+
+    if (BDS_ComplexPacking(bds)) {
+	fprintf(stderr,"*** Cannot decode complex packed fields n=%d***\n", n);
+	exit(8);
+	for (i = 0; i < n; i++) {
+	    *flt++ = UNDEFINED;
+	}
+	return;
+    }
+
+    if (BDS_Harmonic(bds)) {
+        bits = bds + 15;
+        /* fill in global mean */
+        *flt++ = BDS_Harmonic_RefValue(bds);
+        n -= 1; 
+    }
+    else {
+        bits = bds + 11;  
+    }
+
+    tbits = bbits = 0;
+
+    /* assume integer has 32+ bits */
+    if (n_bits <= 25) {
+        jmask = (1 << n_bits) - 1;
+        t_bits = 0;
+
+        if (bitmap) {
+	    for (i = 0; i < n; i++) {
+		/* check bitmap */
+		mask_idx = i & 7;
+		if (mask_idx == 0) bbits = *bitmap++;
+	        if ((bbits & map_masks[mask_idx]) == 0) {
+		    *flt++ = UNDEFINED;
+		    continue;
+	        }
+
+	        while (t_bits < n_bits) {
+	            tbits = (tbits * 256) + *bits++;
+	            t_bits += 8;
+	        }
+	        t_bits -= n_bits;
+	        j = (tbits >> t_bits) & jmask;
+	        *flt++ = ref + scale*j;
+            }
+        }
+        else {
+	    for (i = 0; i < n; i++) {
+		if (n_bits - t_bits > 8) {
+                    tbits = (tbits << 16) | (bits[0] << 8) | (bits[1]);
+		    bits += 2;
+                    t_bits += 16;
+		}
+                while (t_bits < n_bits) {
+                    tbits = (tbits * 256) + *bits++;
+                    t_bits += 8;
+                }
+                t_bits -= n_bits;
+                flt[i] = (tbits >> t_bits) & jmask;
+            }
+	    /* at least this vectorizes :) */
+	    for (i = 0; i < n; i++) {
+		flt[i] = ref + scale*flt[i];
+	    }
+        }
+    }
+    else {
+	/* older unoptimized code, not often used */
+        c_bits = 8;
+        map_mask = 128;
+        while (n-- > 0) {
+	    if (bitmap) {
+	        j = (*bitmap & map_mask);
+	        if ((map_mask >>= 1) == 0) {
+		    map_mask = 128;
+		    bitmap++;
+	        }
+	        if (j == 0) {
+		    *flt++ = UNDEFINED;
+		    continue;
+	        }
+	    }
+
+	    jj = 0.0;
+	    j_bits = n_bits;
+	    while (c_bits <= j_bits) {
+	        if (c_bits == 8) {
+		    jj = jj * 256.0  + (double) (*bits++);
+		    j_bits -= 8;
+	        }
+	        else {
+		    jj = (jj * shift[c_bits]) + (double) (*bits & mask[c_bits]);
+		    bits++;
+		    j_bits -= c_bits;
+		    c_bits = 8;
+	        }
+	    }
+	    if (j_bits) {
+	        c_bits -= j_bits;
+	        jj = (jj * shift[j_bits]) + (double) ((*bits >> c_bits) & mask[j_bits]);
+	    }
+	    *flt++ = ref + scale*jj;
+        }
+    }
+    return;
+}

wgrib/CMakeLists.txt

@@ -5,7 +5,21 @@
 
 add_compile_definitions("DEF_T62_NCEP_TABLE=opn" "FAST_GRIBTAB" "P_TABLE_FIRST")
 
-set(c_src wgrib.c)
+set(c_src BDS_NValues.c	  ectable_132.c  ectable_210.c	nceptab_130.c
+BDSunpk.c	  ectable_133.c  ectable_211.c	nceptab_131.c
+cnames.c	  ectable_140.c  ectable_228.c	nceptab_133.c
+cptectable_254.c  ectable_150.c  ensemble.c	nceptab_140.c
+dwdtable_002.c	  ectable_151.c  flt2ieee.c	nceptab_141.c
+dwdtable_201.c	  ectable_160.c  gds_grid.c	nceptable_opn.c
+dwdtable_202.c	  ectable_162.c  gribtable.c	nceptable_reanal.c
+dwdtable_203.c	  ectable_170.c  ibm2flt.c	PDS_date.c
+dwdtable_204.c	  ectable_171.c  intpower.c	PDStimes.c
+dwdtable_205.c	  ectable_172.c  jra55_200.c	readgrib.c
+ec_ext.c	  ectable_173.c  levels.c	seekgrib.c
+ectable_128.c	  ectable_174.c  mdl_nceptab.c	wgrib_main.c
+ectable_129.c	  ectable_180.c  missing.c	wrtieee.c
+ectable_130.c	  ectable_190.c  nceptab_128.c
+ectable_131.c	  ectable_200.c  nceptab_129.c)
 
 set(exe_name wgrib)
 add_executable(${exe_name} ${c_src})

wgrib/Changes

@@ -1,19 +1,20 @@
-Upgrade from v1.8.0.12g (8/06) to v1.8.0.12o (5/07)
+May 2026
 
-Changes:
+   change code from wgrib.c to source code (*.c *.h)
+   update so that there are no warning with oneapi icx.
+   Note: this is wgrib v1.8.5
 
-        updated ncep grib table 129, 130 to on388
-        updated levels to on388
-        changed units ncep grib table 128 WTPC(186)  K->C
-        support for updates in the grib table by process owners (compile change)
-        increased stdin line buffer to 2000 characters
-        level 210 is now WMO standard for non-NCEP files
+   History: wgrib was written when Unix was a strange operating system.  
+      The various Unix systems had their own versions of "make", and 
+      other operating systems often lacked a make program.  The C 
+      programming language was the new kid on the block, and you learned 
+      from the K&R book.
 
-Reason for upgrade:
+      To help users who didn't understand make and C compilers, the script
+      src2all was written to combine all the *.c and *.h files into wgrib.c
+      so that compiling wgrib would not have to use make.  So compiling wgrib
+      was as simple as:  cc -o wgrib wgrib.c
 
-	Request from marine branch: need updated ncep tables
-
-Comments:
-        No major changes since last update except needed table updates.
-        Compile options change so that table changes can be updated with new JIF.
-	Little change in the run time is expected. 
+      When wgrib was added to the wgrib2 github repository, "wgrib.c" was
+      added.  This makes modifying the code a PITA.  
+