added support for GRIB Gaussian grid parameter N

src/grid.c

@@ -2569,6 +2569,8 @@ int gridGenerate(grid_t grid)
 	if ( grid.xsize > 0 ) gridDefXsize(gridID, grid.xsize);
 	if ( grid.ysize > 0 ) gridDefYsize(gridID, grid.ysize);
 
+        if ( grid.type == GRID_GAUSSIAN ) gridDefNP(gridID, grid.np);
+
 	if ( grid.nvertex > 0 )
 	  gridDefNvertex(gridID, grid.nvertex);
 
@@ -2647,6 +2649,7 @@ int gridGenerate(grid_t grid)
       }
     case GRID_GAUSSIAN_REDUCED:
       {
+	gridDefNP(gridID, grid.np);
 	gridDefYsize(gridID, grid.ysize);
 	gridDefRowlon(gridID, grid.ysize, grid.rowlon);
 

src/stream_cgribex.c

@@ -185,6 +185,7 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid)
 	grid->size  = ISEC4_NumValues;
 	grid->xsize = ISEC2_NumLon;
 	grid->ysize = ISEC2_NumLat;
+        if ( gridtype == GRID_GAUSSIAN ) grid->np = ISEC2_NumPar;
 	grid->xinc  = 0;
 	grid->yinc  = 0;
 	grid->xdef  = 0;
@@ -237,6 +238,7 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid)
       }
     case GRID_GAUSSIAN_REDUCED:
       {
+        grid->np     = ISEC2_NumPar;
 	grid->size   = ISEC4_NumValues;
         grid->rowlon = ISEC2_RowLonPtr;
 	grid->ysize  = ISEC2_NumLat;
@@ -1778,7 +1780,11 @@ void cgribexDefGrid(int *isec1, int *isec2, int *isec4, int gridID)
 	  ISEC2_LonIncr = 0;
 
 	if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
-	  ISEC2_NumPar = nlat/2;
+          {
+            int np = gridInqNP(gridID);
+            if ( np == 0 ) np = nlat/2;
+            ISEC2_NumPar = np;
+          }
 	else
 	  {
 	    ISEC2_LatIncr = NINT(yinc*1000);

src/stream_gribapi.c

@@ -270,9 +270,16 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
 	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
 	nlat = lpar;
 
+	if ( gridtype == GRID_GAUSSIAN )
+          {
+            GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
+            grid->np = lpar;
+          }
+
 	if ( numberOfPoints != nlon*nlat )
 	  Error("numberOfPoints (%d) and gridSize (%d) differ!",
 		(int)numberOfPoints, nlon*nlat);
+
 	grid->size  = numberOfPoints;
 	grid->xsize = nlon;
 	grid->ysize = nlat;
@@ -331,6 +338,9 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
 	size_t dummy;
 	long *pl;
 
+        GRIB_CHECK(grib_get_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar), 0);
+        grid->np = lpar;
+
 	GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
 	nlat = lpar;
 
@@ -2147,7 +2157,11 @@ void gribapiDefGrid(grib_handle *gh, int gridID, int ljpeg)
 	  ISEC2_LonIncr = 0;
 	*/
 	if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
-	  GRIB_CHECK(grib_set_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", nlat/2), 0);
+          {
+            int np = gridInqNP(gridID);
+            if ( np == 0 ) np = nlat/2;
+            GRIB_CHECK(grib_set_long(gh, "numberOfParallelsBetweenAPoleAndTheEquator", np), 0);
+          }
 	else
 	  {
 	    latIncr = yinc;