#include "grib_int.h"
#include <string.h>
#include <math.h>
#include "dmemory.h"

void
gribPrintSec0(int *isec0)
{
  /*

    Print the information in the Indicator
    Section (Section 0) of decoded GRIB data.

    Input Parameters:

       isec0 - Array of decoded integers from Section 0


    Converted from EMOS routine GRPRS0.

       Uwe Schulzweida   MPIfM   01/04/2001

  */

  grsdef();

  fprintf(grprsm, " \n");
  fprintf(grprsm, " Section 0 - Indicator Section.       \n");
  fprintf(grprsm, " -------------------------------------\n");
  fprintf(grprsm, " Length of GRIB message (octets).     %9d\n", ISEC0_GRIB_Len);
  fprintf(grprsm, " GRIB Edition Number.                 %9d\n", ISEC0_GRIB_Version);
}

void
gribPrintSec1(int *isec0, int *isec1)
{
  /*

    Print the information in the Product Definition
    Section (Section 1) of decoded GRIB data.

    Input Parameters
       isec0 - Array of decoded integers from Section 0

       isec1 - Array of decoded integers from Section 1

    Comments:

       When decoding data from Experimental Edition or Edition 0,
       routine GRIBEX adds the additional fields available in
       Edition 1.


    Converted from EMOS routine GRPRS1.

       Uwe Schulzweida   MPIfM   01/04/2001

  */

  int iprev, icurr, ioffset;
  int ibit, ierr, iout, iyear;
  int jiloop;
  float value;
  char hversion[9];

  grsdef();

  /*
    -----------------------------------------------------------------
    Section 0 . Print required information.
    -----------------------------------------------------------------
  */

  fprintf(grprsm, " \n");
  fprintf(grprsm, " Section 1 - Product Definition Section.\n");
  fprintf(grprsm, " ---------------------------------------\n");

  fprintf(grprsm, " Code Table 2 Version Number.         %9d\n", isec1[0]);
  fprintf(grprsm, " Originating centre identifier.       %9d\n", isec1[1]);
  fprintf(grprsm, " Model identification.                %9d\n", isec1[2]);
  fprintf(grprsm, " Grid definition.                     %9d\n", isec1[3]);

  ibit = 8;
  prtbin(isec1[4], ibit, &iout, &ierr);
  fprintf(grprsm, " Flag (Code Table 1)                   %8.8d\n", iout);
  fprintf(grprsm, " Parameter identifier (Code Table 2). %9d\n", isec1[5]);

  if (isec1[5] != 127)
    {
      fprintf(grprsm, " Type of level (Code Table 3).        %9d\n", isec1[6]);
      fprintf(grprsm, " Value 1 of level (Code Table 3).     %9d\n", isec1[7]);
      fprintf(grprsm, " Value 2 of level (Code Table 3).     %9d\n", isec1[8]);
    }
  else
    {
      fprintf(grprsm, " Satellite identifier.                %9d\n", isec1[6]);
      fprintf(grprsm, " Spectral band.                       %9d\n", isec1[7]);
    }

  iyear = isec1[9];
  if (iyear != 255)
    {
      int date, time;
      gribDateTime(isec1, &date, &time);
      iyear = date / 10000;
      fprintf(grprsm, " Year of reference time of data.      %9d  (%4d)\n", isec1[9], iyear);
    }
  else
    {
      fprintf(grprsm, " Year of reference time of data MISSING  (=255)\n");
    }

  fprintf(grprsm, " Month of reference time of data.     %9d\n", isec1[10]);
  fprintf(grprsm, " Day of reference time of data.       %9d\n", isec1[11]);
  fprintf(grprsm, " Hour of reference time of data.      %9d\n", isec1[12]);
  fprintf(grprsm, " Minute of reference time of data.    %9d\n", isec1[13]);
  fprintf(grprsm, " Time unit (Code Table 4).            %9d\n", isec1[14]);
  fprintf(grprsm, " Time range one.                      %9d\n", isec1[15]);
  fprintf(grprsm, " Time range two.                      %9d\n", isec1[16]);
  fprintf(grprsm, " Time range indicator (Code Table 5)  %9d\n", isec1[17]);
  fprintf(grprsm, " Number averaged.                     %9d\n", isec1[18]);
  fprintf(grprsm, " Number missing from average.         %9d\n", isec1[19]);
  /*
     All ECMWF data in GRIB Editions before Edition 1 is decoded
     as 20th century data. Other centres are decoded as missing.
  */
  if (isec0[1] < 1 && isec1[1] != 98)
    fprintf(grprsm, " Century of reference time of data.   Not given\n");
  else
    fprintf(grprsm, " Century of reference time of data.   %9d\n", isec1[20]);

  //   Print sub-centre
  fprintf(grprsm, " Sub-centre identifier.               %9d\n", ISEC1_SubCenterID);

  //   Decimal scale factor
  fprintf(grprsm, " Units decimal scaling factor.        %9d\n", isec1[22]);

  /*
    -----------------------------------------------------------------
    Section 1 . Print local DWD information.
    -----------------------------------------------------------------
  */
  if ((ISEC1_CenterID == 78 || ISEC1_CenterID == 215 || ISEC1_CenterID == 250) && (isec1[36] == 253 || isec1[36] == 254))
    {
      fprintf(grprsm, " DWD local usage identifier.          %9d\n", isec1[36]);
      if (isec1[36] == 253) fprintf(grprsm, " (Database labelling and ensemble forecast)\n");
      if (isec1[36] == 254) fprintf(grprsm, " (Database labelling)\n");

      fprintf(grprsm, " Year of database entry                     %3d  (%4d)\n", isec1[43], 1900 + isec1[43]);
      fprintf(grprsm, " Month of database entry                    %3d\n", isec1[44]);
      fprintf(grprsm, " Day of database entry                      %3d\n", isec1[45]);
      fprintf(grprsm, " Hour of database entry                     %3d\n", isec1[46]);
      fprintf(grprsm, " Minute of database entry                   %3d\n", isec1[47]);
      fprintf(grprsm, " DWD experiment number                %9d\n", isec1[48]);
      fprintf(grprsm, " DWD run type                         %9d\n", isec1[49]);
      if (isec1[36] == 253)
        {
          fprintf(grprsm, " User id                              %9d\n", isec1[50]);
          fprintf(grprsm, " Experiment identifier                %9d\n", isec1[51]);
          fprintf(grprsm, " Ensemble identification type         %9d\n", isec1[52]);
          fprintf(grprsm, " Number of ensemble members           %9d\n", isec1[53]);
          fprintf(grprsm, " Actual number of ensemble member     %9d\n", isec1[54]);
          fprintf(grprsm, " Model version                            %2d.%2.2d\n", isec1[55], isec1[56]);
        }
    }

  /*
    -----------------------------------------------------------------
    Section 2 . Print local ECMWF information.
    -----------------------------------------------------------------
  */
  /*
    Regular MARS labelling, or reformatted Washington EPS products.
  */
  if ((ISEC1_CenterID == 98 && ISEC1_LocalFLag == 1) || (ISEC1_SubCenterID == 98 && ISEC1_LocalFLag == 1)
      || (ISEC1_CenterID == 7 && ISEC1_SubCenterID == 98))
    {
      /*   Parameters common to all definitions.  */

      fprintf(grprsm, " ECMWF local usage identifier.        %9d\n", isec1[36]);
      if (isec1[36] == 1) fprintf(grprsm, " (Mars labelling or ensemble forecast)\n");
      if (isec1[36] == 2) fprintf(grprsm, " (Cluster means and standard deviations)\n");
      if (isec1[36] == 3) fprintf(grprsm, " (Satellite image data)\n");
      if (isec1[36] == 4) fprintf(grprsm, " (Ocean model data)\n");
      if (isec1[36] == 5) fprintf(grprsm, " (Forecast probability data)\n");
      if (isec1[36] == 6) fprintf(grprsm, " (Surface temperature data)\n");
      if (isec1[36] == 7) fprintf(grprsm, " (Sensitivity data)\n");
      if (isec1[36] == 8) fprintf(grprsm, " (ECMWF re-analysis data)\n");
      if (isec1[36] == 9) fprintf(grprsm, " (Singular vectors and ensemble perturbations)\n");
      if (isec1[36] == 10) fprintf(grprsm, " (EPS tubes)\n");
      if (isec1[36] == 11) fprintf(grprsm, " (Supplementary data used by analysis)\n");
      if (isec1[36] == 13) fprintf(grprsm, " (Wave 2D spectra direction and frequency)\n");

      fprintf(grprsm, " Class.                               %9d\n", isec1[37]);
      fprintf(grprsm, " Type.                                %9d\n", isec1[38]);
      fprintf(grprsm, " Stream.                              %9d\n", isec1[39]);
      snprintf(hversion, sizeof(hversion), "%4s", (char *) &isec1[40]);
      hversion[4] = 0;
      fprintf(grprsm, " Version number or Experiment identifier.  %4s\n", hversion);
      /*
        ECMWF Local definition 1.
        (MARS labelling or ensemble forecast data)
      */
      if (isec1[36] == 1)
        {
          fprintf(grprsm, " Forecast number.                     %9d\n", isec1[41]);
          if (isec1[39] != 1090) fprintf(grprsm, " Total number of forecasts.           %9d\n", isec1[42]);

          return;
        }
      /*
        ECMWF Local definition 2.
        (Cluster means and standard deviations)
      */
      if (isec1[36] == 2)
        {
          fprintf(grprsm, " Cluster number.                      %9d\n", isec1[41]);
          fprintf(grprsm, " Total number of clusters.            %9d\n", isec1[42]);
          fprintf(grprsm, " Clustering method.                   %9d\n", isec1[43]);
          fprintf(grprsm, " Start time step when clustering.     %9d\n", isec1[44]);
          fprintf(grprsm, " End time step when clustering.       %9d\n", isec1[45]);
          fprintf(grprsm, " Northern latitude of domain.         %9d\n", isec1[46]);
          fprintf(grprsm, " Western longitude of domain.         %9d\n", isec1[47]);
          fprintf(grprsm, " Southern latitude of domain.         %9d\n", isec1[48]);
          fprintf(grprsm, " Eastern longitude of domain.         %9d\n", isec1[49]);
          fprintf(grprsm, " Operational forecast in cluster      %9d\n", isec1[50]);
          fprintf(grprsm, " Control forecast in cluster          %9d\n", isec1[51]);
          fprintf(grprsm, " Number of forecasts in cluster.      %9d\n", isec1[52]);

          for (int jloop = 0; jloop < isec1[52]; jloop++)
            fprintf(grprsm, " Forecast number                      %9d\n", isec1[jloop + 53]);

          return;
        }
      /*
        ECMWF Local definition 3.
        (Satellite image data)
      */
      if (isec1[36] == 3)
        {
          fprintf(grprsm, " Satellite spectral band.             %9d\n", isec1[41]);
          fprintf(grprsm, " Function code.                       %9d\n", isec1[42]);
          return;
        }
      /*
        ECMWF Local definition 4.
        (Ocean model data)
      */
      if (isec1[36] == 4)
        {
          fprintf(grprsm, " Satellite spectral band.             %9d\n", isec1[41]);
          if (isec1[39] != 1090) fprintf(grprsm, " Function code.                       %9d\n", isec1[42]);
          fprintf(grprsm, " Coordinate structure definition.\n");
          fprintf(grprsm, " Fundamental spatial reference system.%9d\n", isec1[43]);
          fprintf(grprsm, " Fundamental time reference.          %9d\n", isec1[44]);
          fprintf(grprsm, " Space unit flag.                     %9d\n", isec1[45]);
          fprintf(grprsm, " Vertical coordinate definition.      %9d\n", isec1[46]);
          fprintf(grprsm, " Horizontal coordinate definition.    %9d\n", isec1[47]);
          fprintf(grprsm, " Time unit flag.                      %9d\n", isec1[48]);
          fprintf(grprsm, " Time coordinate definition.          %9d\n", isec1[49]);
          fprintf(grprsm, " Position definition.     \n");
          fprintf(grprsm, " Mixed coordinate field flag.         %9d\n", isec1[50]);
          fprintf(grprsm, " Coordinate 1 flag.                   %9d\n", isec1[51]);
          fprintf(grprsm, " Averaging flag.                      %9d\n", isec1[52]);
          fprintf(grprsm, " Position of level 1.                 %9d\n", isec1[53]);
          fprintf(grprsm, " Position of level 2.                 %9d\n", isec1[54]);
          fprintf(grprsm, " Coordinate 2 flag.                   %9d\n", isec1[55]);
          fprintf(grprsm, " Averaging flag.                      %9d\n", isec1[56]);
          fprintf(grprsm, " Position of level 1.                 %9d\n", isec1[57]);
          fprintf(grprsm, " Position of level 2.                 %9d\n", isec1[58]);
          fprintf(grprsm, " Grid Definition.\n");
          fprintf(grprsm, " Coordinate 3 flag (x-axis)           %9d\n", isec1[59]);
          fprintf(grprsm, " Coordinate 4 flag (y-axis)           %9d\n", isec1[60]);
          fprintf(grprsm, " Coordinate 4 of first grid point.    %9d\n", isec1[61]);
          fprintf(grprsm, " Coordinate 3 of first grid point.    %9d\n", isec1[62]);
          fprintf(grprsm, " Coordinate 4 of last grid point.     %9d\n", isec1[63]);
          fprintf(grprsm, " Coordinate 3 of last grid point.     %9d\n", isec1[64]);
          fprintf(grprsm, " i - increment.                       %9d\n", isec1[65]);
          fprintf(grprsm, " j - increment.                       %9d\n", isec1[66]);
          fprintf(grprsm, " Flag for irregular grid coordinates. %9d\n", isec1[67]);
          fprintf(grprsm, " Flag for normal or staggered grids.  %9d\n", isec1[68]);
          fprintf(grprsm, " Further information.\n");
          fprintf(grprsm, " Further information flag.            %9d\n", isec1[69]);
          fprintf(grprsm, " Auxiliary information.\n");
          fprintf(grprsm, " No. entries in horizontal coordinate %9d\n", isec1[70]);
          fprintf(grprsm, " No. entries in mixed coordinate defn.%9d\n", isec1[71]);
          fprintf(grprsm, " No. entries in grid coordinate list. %9d\n", isec1[72]);
          fprintf(grprsm, " No. entries in auxiliary array.      %9d\n", isec1[73]);
          /*
            Horizontal coordinate supplement.
          */
          fprintf(grprsm, " Horizontal coordinate supplement.\n");
          if (isec1[70] == 0)
            {
              fprintf(grprsm, "(None).\n");
            }
          else
            {
              fprintf(grprsm, "Number of items = %d\n", isec1[70]);
              for (int jloop = 0; jloop < isec1[70]; jloop++) fprintf(grprsm, "         %12d\n", isec1[74 + jloop]);
            }
          /*
            Mixed coordinate definition.
          */
          fprintf(grprsm, " Mixed coordinate definition.\n");
          if (isec1[71] == 0)
            {
              fprintf(grprsm, "(None).\n");
            }
          else
            {
              fprintf(grprsm, "Number of items = %d\n", isec1[71]);
              ioffset = 74 + isec1[70];
              for (int jloop = 0; jloop < isec1[71]; jloop++) fprintf(grprsm, "         %12d\n", isec1[ioffset + jloop]);
            }
          /*
            Grid coordinate list.
          */
          fprintf(grprsm, " Grid coordinate list. \n");
          if (isec1[72] == 0)
            {
              fprintf(grprsm, "(None).\n");
            }
          else
            {
              fprintf(grprsm, "Number of items = %d\n", isec1[72]);
              ioffset = 74 + isec1[70] + isec1[71];
              for (int jloop = 0; jloop < isec1[72]; jloop++) fprintf(grprsm, "         %12d\n", isec1[ioffset + jloop]);
            }
          /*
            Auxiliary array.
          */
          fprintf(grprsm, " Auxiliary array.      \n");
          if (isec1[73] == 0)
            {
              fprintf(grprsm, "(None).\n");
            }
          else
            {
              fprintf(grprsm, "Number of items = %d\n", isec1[73]);
              ioffset = 74 + isec1[70] + isec1[71] + isec1[72];
              for (int jloop = 0; jloop < isec1[73]; jloop++) fprintf(grprsm, "         %12d\n", isec1[ioffset + jloop]);
            }
          /*
            Post-auxiliary array.
          */
          fprintf(grprsm, " Post-auxiliary array. \n");
          ioffset = 74 + isec1[70] + isec1[71] + isec1[72] + isec1[73];
          if (isec1[ioffset] == 0)
            {
              fprintf(grprsm, "(None).\n");
            }
          else
            {
              fprintf(grprsm, "Number of items = %d\n", isec1[ioffset]);
              for (int jloop = 1; jloop < isec1[ioffset]; jloop++) fprintf(grprsm, "         %12d\n", isec1[ioffset + jloop]);
            }

          return;
        }
      /*
        ECMWF Local definition 5.
        (Forecast probability data)
      */
      if (isec1[36] == 5)
        {
          fprintf(grprsm, " Forecast probability number          %9d\n", isec1[41]);
          fprintf(grprsm, " Total number of forecast probabilities %7d\n", isec1[42]);
          fprintf(grprsm, " Threshold units decimal scale factor %9d\n", isec1[43]);
          fprintf(grprsm, " Threshold indicator(1=lower,2=upper,3=both) %2d\n", isec1[44]);
          if (isec1[44] != 2) fprintf(grprsm, " Lower threshold value                %9d\n", isec1[45]);
          if (isec1[44] != 1) fprintf(grprsm, " Upper threshold value                %9d\n", isec1[46]);
          return;
        }
      /*
        ECMWF Local definition 6.
        (Surface temperature data)
      */
      if (isec1[36] == 6)
        {
          iyear = isec1[43];
          if (iyear > 100)
            {
              if (iyear < 19000000) iyear = iyear + 19000000;
              fprintf(grprsm, " Date of SST field used               %9d\n", iyear);
            }
          else
            fprintf(grprsm, "Date of SST field used               Not given\n");
        }
      if (isec1[44] == 0) fprintf(grprsm, " Type of SST field (= climatology)    %9d\n", isec1[44]);
      if (isec1[44] == 1) fprintf(grprsm, " Type of SST field (= 1/1 degree)     %9d\n", isec1[44]);
      if (isec1[44] == 2) fprintf(grprsm, " Type of SST field (= 2/2 degree)     %9d\n", isec1[44]);

      fprintf(grprsm, " Number of ICE fields used:           %9d\n", isec1[45]);

      for (int jloop = 1; jloop <= isec1[45]; jloop++)
        {
          iyear = isec1[44 + (jloop * 2)];
          if (iyear > 100)
            {
              if (iyear < 19000000) iyear = iyear + 19000000;
              fprintf(grprsm, " Date of ICE field%3d                 %9d\n", jloop, iyear);
              fprintf(grprsm, " Satellite number (ICE field%3d)      %9d\n", jloop, isec1[45 + (jloop * 2)]);
            }
          else
            fprintf(grprsm, "Date of SST field used               Not given\n");
        }
      /*
        ECMWF Local definition 7.
        (Sensitivity data)
      */
      if (isec1[36] == 7)
        {
          if (isec1[38] == 51) fprintf(grprsm, " Forecast number                      %9d\n", isec1[41]);
          if (isec1[38] != 51) fprintf(grprsm, " Iteration number                     %9d\n", isec1[41]);
          if (isec1[38] != 52) fprintf(grprsm, " Total number of diagnostics          %9d\n", isec1[42]);
          if (isec1[38] == 52) fprintf(grprsm, " No.interations in diag. minimisation %9d\n", isec1[42]);
          fprintf(grprsm, " Domain(0=Global,1=Europe,2=N.Hem.,3=S.Hem.) %2d\n", isec1[43]);
          fprintf(grprsm, " Diagnostic number                    %9d\n", isec1[44]);
        }
      /*
        ECMWF Local definition 8.
        (ECMWF re-analysis data)
      */
      if (isec1[36] == 8)
        {
          if ((isec1[39] == 1043) || (isec1[39] == 1070) || (isec1[39] == 1071))
            {
              fprintf(grprsm, " Interval between reference times     %9d\n", isec1[41]);
              for (int jloop = 43; jloop <= 54; jloop++)
                {
                  jiloop = jloop + 8;
                  fprintf(grprsm, " ERA section 1 octet %2d.              %9d\n", jiloop, isec1[jloop - 1]);
                }
            }
          else
            {
              for (int jloop = 42; jloop <= 54; jloop++)
                {
                  jiloop = jloop + 8;
                  fprintf(grprsm, " ERA section 1 octet %2d.              %9d\n", jiloop, isec1[jloop - 1]);
                }
            }
          return;
        }

      if (isec1[38] > 4 && isec1[38] < 9)
        {
          fprintf(grprsm, " Simulation number.                   %9d\n", isec1[41]);
          fprintf(grprsm, " Total number of simulations.         %9d\n", isec1[42]);
        }
      /*
        ECMWF Local definition 9.
        (Singular vectors and ensemble perturbations)
      */
      if (isec1[36] == 9)
        {
          if (isec1[38] == 60) fprintf(grprsm, " Perturbed ensemble forecast number   %9d\n", isec1[41]);
          if (isec1[38] == 61) fprintf(grprsm, " Initial state perturbation number    %9d\n", isec1[41]);
          if (isec1[38] == 62) fprintf(grprsm, " Singular vector number               %9d\n", isec1[41]);
          if (isec1[38] == 62)
            {
              fprintf(grprsm, " Number of iterations                 %9d\n", isec1[42]);
              fprintf(grprsm, " Number of singular vectors computed  %9d\n", isec1[43]);
              fprintf(grprsm, " Norm used at initial time            %9d\n", isec1[44]);
              fprintf(grprsm, " Norm used at final time              %9d\n", isec1[45]);
              fprintf(grprsm, " Multiplication factor                %9d\n", isec1[46]);
              fprintf(grprsm, " Latitude of north-west corner        %9d\n", isec1[47]);
              fprintf(grprsm, " Longitude of north-west corner       %9d\n", isec1[48]);
              fprintf(grprsm, " Latitude of south-east corner        %9d\n", isec1[49]);
              fprintf(grprsm, " Longitude of south-east corner       %9d\n", isec1[50]);
              fprintf(grprsm, " Accuracy                             %9d\n", isec1[51]);
              fprintf(grprsm, " Number of singular vectors evolved   %9d\n", isec1[52]);
              fprintf(grprsm, " Ritz number one                      %9d\n", isec1[53]);
              fprintf(grprsm, " Ritz number two                      %9d\n", isec1[54]);
            }
        }
      /*
        ECMWF Local definition 10.
        (EPS tubes)
      */
      if (isec1[36] == 10)
        {
          fprintf(grprsm, " Tube number                          %9d\n", isec1[41]);
          fprintf(grprsm, " Total number of tubes                %9d\n", isec1[42]);
          fprintf(grprsm, " Central cluster definition           %9d\n", isec1[43]);
          fprintf(grprsm, " Parameter                            %9d\n", isec1[44]);
          fprintf(grprsm, " Type of level                        %9d\n", isec1[45]);
          fprintf(grprsm, " Northern latitude of domain of tubing%9d\n", isec1[46]);
          fprintf(grprsm, " Western longitude of domain of tubing%9d\n", isec1[47]);
          fprintf(grprsm, " Southern latitude of domain of tubing%9d\n", isec1[48]);
          fprintf(grprsm, " Eastern longitude of domain of tubing%9d\n", isec1[49]);
          fprintf(grprsm, " Tube number of operational forecast  %9d\n", isec1[50]);
          fprintf(grprsm, " Tube number of control forecast      %9d\n", isec1[51]);
          fprintf(grprsm, " Height/pressure of level             %9d\n", isec1[52]);
          fprintf(grprsm, " Reference step                       %9d\n", isec1[53]);
          fprintf(grprsm, " Radius of central cluster            %9d\n", isec1[54]);
          fprintf(grprsm, " Ensemble standard deviation          %9d\n", isec1[55]);
          fprintf(grprsm, " Dist.of tube extreme to ensemble mean%9d\n", isec1[56]);
          fprintf(grprsm, " Number of forecasts in the tube      %9d\n", isec1[57]);

          fprintf(grprsm, " List of ensemble forecast numbers:\n");
          for (int jloop = 1; jloop <= isec1[57]; jloop++) fprintf(grprsm, "    %9d\n", isec1[57 + jloop]);
        }
      /*
        ECMWF Local definition 11.
        (Supplementary data used by the analysis)
      */
      if (isec1[36] == 11)
        {
          fprintf(grprsm, " Details of analysis which used the supplementary data:\n");
          fprintf(grprsm, "   Class                              %9d\n", isec1[41]);
          fprintf(grprsm, "   Type                               %9d\n", isec1[42]);
          fprintf(grprsm, "   Stream                             %9d\n", isec1[43]);
          /*
          snprintf(hversion, sizeof(hversion), "%8d", isec1[44]);
          fprintf(grprsm, "   Version number/experiment identifier:   %4s\n", &hversion[4]);
          */
          iyear = isec1[45];
          iyear = iyear + ((iyear > 50) ? 1900 : 2000);

          fprintf(grprsm, "   Year                               %9d\n", iyear);
          fprintf(grprsm, "   Month                              %9d\n", isec1[46]);
          fprintf(grprsm, "   Day                                %9d\n", isec1[47]);
          fprintf(grprsm, "   Hour                               %9d\n", isec1[48]);
          fprintf(grprsm, "   Minute                             %9d\n", isec1[49]);
          fprintf(grprsm, "   Century                            %9d\n", isec1[50]);
          fprintf(grprsm, "   Originating centre                 %9d\n", isec1[51]);
          fprintf(grprsm, "   Sub-centre                         %9d\n", isec1[52]);
        }
      /*
        ECMWF Local definition 12.
      */
      if (isec1[36] == 12)
        {
          fprintf(grprsm, " (Mean, average, etc)\n");
          fprintf(grprsm, " Start date of the period              %8d\n", isec1[41]);
          fprintf(grprsm, " Start time of the period                  %4.4d\n", isec1[42]);
          fprintf(grprsm, " Finish date of the period             %8d\n", isec1[43]);
          fprintf(grprsm, " Finish time of the period                 %4.4d\n", isec1[44]);
          fprintf(grprsm, " Verifying date of the period          %8d\n", isec1[45]);
          fprintf(grprsm, " Verifying time of the period              %4.4d\n", isec1[46]);
          fprintf(grprsm, " Code showing method                   %8d\n", isec1[47]);
          fprintf(grprsm, " Number of different time intervals used  %5d\n", isec1[48]);
          fprintf(grprsm, " List of different time intervals used:\n");
          iprev = isec1[49];
          unsigned icount = 0;
          for (int jloop = 1; jloop <= isec1[48]; jloop++)
            {
              icurr = isec1[48 + jloop];
              if (icurr != iprev)
                {
                  if (icount == 1) fprintf(grprsm, "  - interval %5.4d used       once\n", iprev);
                  if (icount == 2) fprintf(grprsm, "  - interval %5.4d used       twice\n", iprev);
                  if (icount > 2) fprintf(grprsm, "  - interval %5.4d used %5u times\n", iprev, icount);
                  iprev = icurr;
                  icount = 1;
                }
              else
                icount = icount + 1;
            }
          if (icount == 1) fprintf(grprsm, "  - interval %5.4d used       once\n", iprev);
          if (icount == 2) fprintf(grprsm, "  - interval %5.4d used       twice\n", iprev);
          if (icount > 2) fprintf(grprsm, "  - interval %5.4d used %5u times\n", iprev, icount);
        }
      /*
        ECMWF Local definition 13.
        (Wave 2D spectra direction and frequency)
      */
      if (isec1[36] == 13)
        {
          fprintf(grprsm, " Direction number                     %9d\n", isec1[43]);
          fprintf(grprsm, " Frequency number                     %9d\n", isec1[44]);
          fprintf(grprsm, " Total number of directions           %9d\n", isec1[45]);
          fprintf(grprsm, " Total number of frequencies          %9d\n", isec1[46]);
          fprintf(grprsm, " Scale factor applied to directions   %9d\n", isec1[47]);
          fprintf(grprsm, " Scale factor applied to frequencies  %9d\n", isec1[48]);
          fprintf(grprsm, " List of directions:\n");
          for (int jloop = 1; jloop <= isec1[45]; jloop++)
            {
              value = (float) (isec1[48 + jloop]) / (float) (isec1[47]);
              if (isec1[43] == jloop)
                fprintf(grprsm, " %2.2d:%15.7f   <-- this field value\n", jloop, value);
              else
                fprintf(grprsm, "%2.2d:%15.7f\n", jloop, value);
            }
          fprintf(grprsm, " List of frequencies:\n");
          for (int jloop = 1; jloop <= isec1[46]; jloop++)
            {
              value = (float) (isec1[48 + isec1[45] + jloop]) / (float) (isec1[48]);
              if (isec1[44] == jloop)
                fprintf(grprsm, " %2.2d:%15.7f   <-- this field value\n", jloop, value);
              else
                fprintf(grprsm, "%2.2d:%15.7f\n", jloop, value);

              if (isec1[49 + isec1[45] + isec1[46]] != 0)
                {
                  fprintf(grprsm, " System number (65535 = missing)      %9d\n", isec1[49 + isec1[45] + isec1[46]]);
                  fprintf(grprsm, " Method number (65535 = missing)      %9d\n", isec1[50 + isec1[45] + isec1[46]]);
                }
            }
          /*
            ECMWF Local definition 14.
            (Brightness temperature)
          */
          if (isec1[36] == 14)
            {
              fprintf(grprsm, " Channel number                       %9d\n", isec1[43]);
              fprintf(grprsm, " Scale factor applied to frequencies  %9d\n", isec1[44]);
              fprintf(grprsm, " Total number of frequencies          %9d\n", isec1[45]);
              fprintf(grprsm, " List of frequencies:\n");
              for (int jloop = 1; jloop <= isec1[45]; jloop++)
                {
                  value = (float) (isec1[45 + jloop]) / (float) (isec1[44]);
                  if (isec1[43] == jloop)
                    fprintf(grprsm, " %3d:%15.9f   <-- this channel\n", jloop, value);
                  else
                    fprintf(grprsm, " %3d:%15.9f\n", jloop, value);
                }
            }
          /*
            ECMWF Local definition 15.
            (Ocean ensemble seasonal forecast)
          */
          if (isec1[36] == 15)
            {
              fprintf(grprsm, " Ensemble member number               %9d\n", isec1[41]);
              fprintf(grprsm, " System number                        %9d\n", isec1[42]);
              fprintf(grprsm, " Method number                        %9d\n", isec1[43]);
            }
          /*
            ECMWF Local definition 16.
            (Seasonal forecast monthly mean atmosphere data)
          */
          if (isec1[36] == 16)
            {
              fprintf(grprsm, " Ensemble member number               %9d\n", isec1[41]);
              fprintf(grprsm, " System number                        %9d\n", isec1[43]);
              fprintf(grprsm, " Method number                        %9d\n", isec1[44]);
              fprintf(grprsm, " Verifying month                      %9d\n", isec1[45]);
              fprintf(grprsm, " Averaging period                     %9d\n", isec1[46]);
            }
          /*
            ECMWF Local definition 17.
            (Sst or sea-ice used by analysis)
          */
          if (isec1[36] == 17)
            {
              iyear = isec1[43];
              if (iyear > 100)
                {
                  if (iyear < 19000000) iyear = iyear + 19000000;
                  fprintf(grprsm, " Date of sst/ice field used           %9d\n", iyear);
                }
              else
                fprintf(grprsm, " Date of sst/ice field used           Not given\n");

              if (isec1[44] == 0) fprintf(grprsm, " Type of sst/ice field (= climatology)%9d\n", isec1[44]);
              if (isec1[44] == 1) fprintf(grprsm, " Type of sst/ice field (= 1/1 degree) %9d\n", isec1[44]);
              if (isec1[44] == 2) fprintf(grprsm, " Type of sst/ice field (= 2/2 degree) %9d\n", isec1[44]);

              fprintf(grprsm, " Number of ICE fields used:           %9d\n", isec1[45]);

              for (int jloop = 1; jloop < isec1[45]; jloop++)
                {
                  iyear = isec1[44 + (jloop * 2)];
                  if (iyear > 100)
                    {
                      if (iyear < 19000000) iyear = iyear + 19000000;
                      fprintf(grprsm, " Date of ICE field%3d                 %9d\n", jloop, iyear);
                      fprintf(grprsm, " Satellite number (ICE field%3d)      %9d\n", jloop, isec1[45 + (jloop * 2)]);
                    }
                  else
                    fprintf(grprsm, "Date of sst/ice field used           Not given\n");
                }
            }
        }
    }
  /*
    -----------------------------------------------------------------
    Section 3 . Print Washington ensemble product information.
    -----------------------------------------------------------------
  */
  /*
    Washington EPS products (but not reformatted Washington EPS
    products.
  */
  if ((isec1[1] == 7 && isec1[23] == 1) && (!(ISEC1_SubCenterID == 98)))
    {
      /*   CALL KWPRS1 (iSEC0,iSEC1)*/
    }
  /*
    -----------------------------------------------------------------
    Section 4 . Print local MPIM information.
    -----------------------------------------------------------------
  */
  if (isec1[1] == 252 && isec1[36] == 1)
    {
      fprintf(grprsm, " MPIM local usage identifier.         %9d\n", isec1[36]);
      fprintf(grprsm, " Type of ensemble forecast            %9d\n", isec1[37]);
      fprintf(grprsm, " Individual ensemble member           %9d\n", isec1[38]);
      fprintf(grprsm, " Number of forecasts in ensemble      %9d\n", isec1[39]);
    }
}

static void
printQuasi(int *isec2)
{
  /*

    Print the qusai-regular information in the Grid Description
    Section (Section 2) of decoded GRIB data.

    Input Parameters:

       isec2 - Array of decoded integers from Section 2.

    Comments:

       Only data representation types catered for are Gaussian
       grid, latitude/longitude grid, Spherical Harmonics,
       Polar stereographic and Space view perspective.

    Converted from EMOS routine PTQUASI.

       Uwe Schulzweida   MPIfM   01/04/2001

  */

  char yout[64];

  /*
    -----------------------------------------------------------------
    Section 1. Print quasi-grid data.
    -----------------------------------------------------------------
  */
  // See if scanning is north->south or south->north
  fprintf(grprsm, "  Number of points along a parallel varies.\n");

  int ntos = (fmod((double) isec2[10], 128.) < 64);
  fprintf(grprsm, "  Number of points.   Parallel. %s\n", ntos ? "(North to South)" : "(South to North)");

  // Display number of points for each latitude
  int latcnt = isec2[2];
  int nextlat = 0;
  memset(yout, ' ', (size_t) 11);

  for (int j = 0; j < latcnt; ++j)
    {
      nextlat = nextlat + 1;
      snprintf(yout, sizeof(yout), "%4d", nextlat);

      // Finished?
      if (nextlat > latcnt) break;
      if (nextlat == latcnt)
        {
          fprintf(grprsm, " %5d                %-12s\n", isec2[nextlat + 21], yout);
          break;
        }
      // Look for neighbouring latitudes with same number of points
      unsigned nrepeat = 0;

    LABEL110:
      // If neighbouring latitudes have same number of points increase the repeat count.
      if (isec2[nextlat + 21 + 1] == isec2[nextlat + 21])
        {
          nrepeat = nrepeat + 1;
          nextlat = nextlat + 1;
          if (nextlat < latcnt) goto LABEL110;
        }
      // Display neighbouring latitudes with same number of points as 'nn to mm'.
      if (nrepeat >= 1) snprintf(yout + 4, sizeof(yout) - 4, "to %5d", nextlat);
      fprintf(grprsm, " %5d                %-12s\n", isec2[nextlat + 21], yout);
      memset(yout, ' ', (size_t) 11);
    }
}

void
gribPrintSec2DP(int *isec0, int *isec2, double *fsec2)
{
  /*

    Print the information in the Grid Description
    Section (Section 2) of decoded GRIB data.

    Input Parameters:

       isec0  - Array of decoded integers from Section 0

       isec2  - Array of decoded integers from Section 2

       fsec2  - Array of decoded floats from Section 2

    Comments:

       Only data representation types catered for are Gaussian
       grid, latitude/longitude grid, Spherical Harmonics,
       Polar stereographic and Space view perspective.


    Converted from EMOS routine GRPRS2.

       Uwe Schulzweida   MPIfM   01/04/2001

  */

  int ibit, iedit, ierr, iout, iresol;

  grsdef();
  /*
    -----------------------------------------------------------------
    Section 1 . Print GRIB Edition number.
    -----------------------------------------------------------------
  */
  iedit = isec0[1];
  fprintf(grprsm, " \n");
  fprintf(grprsm, " Section 2 - Grid Description Section.\n");
  fprintf(grprsm, " -------------------------------------\n");
  /*
    -----------------------------------------------------------------
    Section 2 . Print spherical harmonic data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 50 || isec2[0] == 60 || isec2[0] == 70 || isec2[0] == 80)
    {
      fprintf(grprsm, " Data represent type = spectral     (Table 6) %9d\n", isec2[0]);
      fprintf(grprsm, " J - Pentagonal resolution parameter.         %9d\n", isec2[1]);
      fprintf(grprsm, " K - Pentagonal resolution parameter.         %9d\n", isec2[2]);
      fprintf(grprsm, " M - Pentagonal resolution parameter.         %9d\n", isec2[3]);
      fprintf(grprsm, " Representation type (Table 9)                %9d\n", isec2[4]);
      fprintf(grprsm, " Representation mode (Table 10).              %9d\n", isec2[5]);
      for (int i = 7; i <= 11; ++i) fprintf(grprsm, " Not used.                                    %9d\n", isec2[i - 1]);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 3 . Print Gaussian grid data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 4 || isec2[0] == 14 || isec2[0] == 24 || isec2[0] == 34)
    {
      fprintf(grprsm, " (Southern latitudes and Western longitudes are negative.)\n");
      fprintf(grprsm, " Data represent type = gaussian     (Table 6) %9d\n", isec2[0]);
      /*
        Quasi-regular grids introduced in Edition 1.
      */
      if (isec2[16] == 0 || iedit < 1)
        fprintf(grprsm, " Number of points along a parallel.           %9d\n", isec2[1]);
      else
        printQuasi(isec2);

      fprintf(grprsm, " Number of points along a meridian.           %9d\n", isec2[2]);
      fprintf(grprsm, " Latitude of first grid point.                %9d\n", isec2[3]);
      fprintf(grprsm, " Longitude of first grid point.               %9d\n", isec2[4]);

      ibit = 8;
      iresol = isec2[5] + isec2[17] + isec2[18];
      prtbin(iresol, ibit, &iout, &ierr);

      fprintf(grprsm, " Resolution and components flag.               %8.8d\n", iout);
      fprintf(grprsm, " Latitude of last grid point.                 %9d\n", isec2[6]);
      fprintf(grprsm, " Longitude of last grid point.                %9d\n", isec2[7]);
      /*
        Print increment if given.
      */
      if (isec2[5] == 128)
        fprintf(grprsm, " i direction (East-West) increment.           %9d\n", isec2[8]);
      else
        fprintf(grprsm, " i direction (East-West) increment            Not given\n");

      fprintf(grprsm, " Number of parallels between pole and equator.%9d\n", isec2[9]);

      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);

      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 4 . Print Latitude / longitude grid data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 0 || isec2[0] == 10 || isec2[0] == 20 || isec2[0] == 30)
    {
      fprintf(grprsm, " (Southern latitudes and Western longitudes are negative.)\n");
      fprintf(grprsm, " Data represent type = lat/long     (Table 6) %9d\n", isec2[0]);
      /*
        Quasi-regular lat/long grids also possible.
      */
      if (isec2[16] == 0)
        fprintf(grprsm, " Number of points along a parallel.           %9d\n", isec2[1]);
      else
        printQuasi(isec2);

      fprintf(grprsm, " Number of points along a meridian.           %9d\n", isec2[2]);
      fprintf(grprsm, " Latitude of first grid point.                %9d\n", isec2[3]);
      fprintf(grprsm, " Longitude of first grid point.               %9d\n", isec2[4]);

      ibit = 8;
      iresol = isec2[5] + isec2[17] + isec2[18];
      prtbin(iresol, ibit, &iout, &ierr);

      fprintf(grprsm, " Resolution and components flag.               %8.8d\n", iout);
      fprintf(grprsm, " Latitude of last grid point.                 %9d\n", isec2[6]);
      fprintf(grprsm, " Longitude of last grid point.                %9d\n", isec2[7]);
      /*
        Print increment if given.
      */
      if (isec2[8] < 0)
        fprintf(grprsm, " i direction (East-West) increment            Not given\n");
      else
        fprintf(grprsm, " i direction (East-West) increment.           %9d\n", isec2[8]);

      if (isec2[9] < 0)
        fprintf(grprsm, " j direction (North-South) increment          Not given\n");
      else
        fprintf(grprsm, " j direction (North-South) increment.         %9d\n", isec2[9]);

      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);

      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 5 . Print polar stereographic data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 5)
    {
      fprintf(grprsm, " (Southern latitudes and Western longitudes are negative.)\n");
      fprintf(grprsm, " Data represent type = polar stereo (Table 6) %9d\n", isec2[0]);
      fprintf(grprsm, " Number of points along X axis.               %9d\n", isec2[1]);
      fprintf(grprsm, " Number of points along Y axis.               %9d\n", isec2[2]);
      fprintf(grprsm, " Latitude of first grid point.                %9d\n", isec2[3]);
      fprintf(grprsm, " Longitude of first grid point.               %9d\n", isec2[4]);
      ibit = 8;
      iresol = isec2[17] + isec2[18];
      prtbin(iresol, ibit, &iout, &ierr);
      fprintf(grprsm, " Resolution and components flag.               %8.8d\n", iout);
      fprintf(grprsm, " Orientation of the grid.                     %9d\n", isec2[6]);
      fprintf(grprsm, " X direction increment.                       %9d\n", isec2[8]);
      fprintf(grprsm, " Y direction increment.                       %9d\n", isec2[9]);
      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);
      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      fprintf(grprsm, " Projection centre flag.                      %9d\n", isec2[12]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 6 . Print Lambert conformal data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 3)
    {
      fprintf(grprsm, " (Southern latitudes and Western longitudes are negative.)\n");
      fprintf(grprsm, " Data represent type = Lambert      (Table 6) %9d\n", isec2[0]);
      fprintf(grprsm, " Number of points along X axis.               %9d\n", isec2[1]);
      fprintf(grprsm, " Number of points along Y axis.               %9d\n", isec2[2]);
      fprintf(grprsm, " Latitude of first grid point.                %9d\n", isec2[3]);
      fprintf(grprsm, " Longitude of first grid point.               %9d\n", isec2[4]);
      ibit = 8;
      iresol = isec2[17] + isec2[18] + isec2[5];
      prtbin(iresol, ibit, &iout, &ierr);
      fprintf(grprsm, " Resolution and components flag.               %8.8d\n", iout);
      fprintf(grprsm, " Orientation of the grid.                     %9d\n", isec2[6]);
      fprintf(grprsm, " X direction increment.                       %9d\n", isec2[8]);
      fprintf(grprsm, " Y direction increment.                       %9d\n", isec2[9]);
      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);
      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      fprintf(grprsm, " Projection centre flag.                      %9d\n", isec2[12]);
      fprintf(grprsm, " Latitude intersection 1 - Latin 1 -.         %9d\n", isec2[13]);
      fprintf(grprsm, " Latitude intersection 2 - Latin 2 -.         %9d\n", isec2[14]);
      fprintf(grprsm, " Latitude of Southern Pole.                   %9d\n", isec2[19]);
      fprintf(grprsm, " Longitude of Southern Pole.                  %9d\n", isec2[20]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 7 . Print space view perspective or orthographic data.
    -----------------------------------------------------------------
  */
  if (isec2[0] == 90)
    {
      fprintf(grprsm, " (Southern latitudes and Western longitudes are negative.)\n");
      fprintf(grprsm, " Data represent type = space/ortho  (Table 6) %9d\n", isec2[0]);
      fprintf(grprsm, " Number of points along X axis.               %9d\n", isec2[1]);
      fprintf(grprsm, " Number of points along Y axis.               %9d\n", isec2[2]);
      fprintf(grprsm, " Latitude of sub-satellite point.             %9d\n", isec2[3]);
      fprintf(grprsm, " Longitude of sub-satellite point.            %9d\n", isec2[4]);
      // iresol = isec2[17] + isec2[18];
      fprintf(grprsm, " Diameter of the earth in x direction.        %9d\n", isec2[6]);
      fprintf(grprsm, " Y coordinate of sub-satellite point.         %9d\n", isec2[9]);
      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);
      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      fprintf(grprsm, " Orientation of the grid.                     %9d\n", isec2[6]);
      fprintf(grprsm, " Altitude of the camera.                      %9d\n", isec2[13]);
      fprintf(grprsm, " Y coordinate of origin of sector image.      %9d\n", isec2[14]);
      fprintf(grprsm, " X coordinate of origin of sector image.      %9d\n", isec2[15]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Section 7.5 . Print ocean data
    -----------------------------------------------------------------
  */
  /*
  if ( isec2[0] == 192 && ISEC1_CenterID == 98 )
    {
      fprintf(grprsm, " Data represent type = ECMWF ocean  (Table 6) %9d\n", isec2[0]);
      if ( isec2[1] ==  32767 )
        fprintf(grprsm, " Number of points along the first axis.       Not used\n");
      else
        fprintf(grprsm, " Number of points along the first axis.       %9d\n", isec2[1]);

      if ( isec2[2] ==  32767 )
        fprintf(grprsm, " Number of points along the second axis.      Not used\n");
      else
        fprintf(grprsm, " Number of points along the second axis.      %9d\n", isec2[2]);

      ibit = 8;
      prtbin(isec2[10], ibit, &iout, &ierr);
      fprintf(grprsm, " Scanning mode flags (Code Table 8)            %8.8d\n", iout);
      goto LABEL800;
    }
    */
  /*
    -----------------------------------------------------------------
    Section 7.6 . Print triangular data
    -----------------------------------------------------------------
  */
  if (isec2[0] == 192 /* && ISEC1_CenterID == 78 */)
    {
      fprintf(grprsm, " Data represent type = triangular   (Table 6) %9d\n", isec2[0]);
      fprintf(grprsm, " Number of factor 2 in factorisation of Ni.   %9d\n", isec2[1]);
      fprintf(grprsm, " Number of factor 3 in factorisation of Ni.   %9d\n", isec2[2]);
      fprintf(grprsm, " Number of diamonds (Nd).                     %9d\n", isec2[3]);
      fprintf(grprsm, " Number of triangular subdivisions of the\n");
      fprintf(grprsm, "           icosahedron (Ni).                  %9d\n", isec2[4]);
      fprintf(grprsm, " Flag for orientation of diamonds (Table A).  %9d\n", isec2[5]);
      fprintf(grprsm, " Latitude of pole point.                      %9d\n", isec2[6]);
      fprintf(grprsm, " Longitude of pole point.                     %9d\n", isec2[7]);
      fprintf(grprsm, " Longitude of the first diamond.              %9d\n", isec2[8]);
      fprintf(grprsm, " Flag for storage sequence (Table B).         %9d\n", isec2[9]);
      fprintf(grprsm, " Number of vertical coordinate parameters.    %9d\n", isec2[11]);
      goto LABEL800;
    }
  /*
    -----------------------------------------------------------------
    Drop through to here => representation type not catered for.
    -----------------------------------------------------------------
  */
  fprintf(grprsm, "GRPRS2 :Data representation type not catered for -%d\n", isec2[0]);

  goto LABEL900;
  /*
    -----------------------------------------------------------------
    Section 8 . Print vertical coordinate parameters,
                rotated grid information,
                stretched grid information, if any.
    -----------------------------------------------------------------
  */
LABEL800:;
  /*
    Vertical coordinate parameters ...
  */
  if (isec2[11] != 0)
    {
      fprintf(grprsm, " \n");
      fprintf(grprsm, " Vertical Coordinate Parameters.\n");
      fprintf(grprsm, " -------------------------------\n");
      for (int i = 10; i < isec2[11] + 10; ++i) fprintf(grprsm, "    %20.12f\n", fsec2[i]);
    }
  /*
    Rotated and stretched grids introduced in Edition 1.
  */
  if (iedit < 1) goto LABEL900;
  /*
    Rotated grid information ...
  */
  if (isec2[0] == 10 || isec2[0] == 30 || isec2[0] == 14 || isec2[0] == 34 || isec2[0] == 60 || isec2[0] == 80 || isec2[0] == 30)
    {
      fprintf(grprsm, " \n");
      fprintf(grprsm, " Latitude of southern pole of rotation.       %9d\n", isec2[12]);
      fprintf(grprsm, " Longitude of southern pole of rotation.      %9d\n", isec2[13]);
      fprintf(grprsm, " Angle of rotation.                     %20.10f\n", fsec2[0]);
    }
  /*
    Stretched grid information ...
  */
  if (isec2[0] == 20 || isec2[0] == 30 || isec2[0] == 24 || isec2[0] == 34 || isec2[0] == 70 || isec2[0] == 80)
    {
      fprintf(grprsm, " \n");
      fprintf(grprsm, " Latitude of pole of stretching.              %9d\n", isec2[14]);
      fprintf(grprsm, " Longitude of pole of stretching.             %9d\n", isec2[15]);
      fprintf(grprsm, " Stretching factor.                     %20.10f\n", fsec2[1]);
    }

LABEL900:;

  return;
}

void
gribPrintSec2SP(int *isec0, int *isec2, float *fsec2sp)
{
  int inum = 10 + isec2[11];

  double *fsec2 = (double *) Malloc((size_t) inum * sizeof(double));
  if (fsec2 == NULL) SysError("No Memory!");

  for (int j = 0; j < inum; ++j) fsec2[j] = fsec2sp[j];

  gribPrintSec2DP(isec0, isec2, fsec2);

  Free(fsec2);
}

void
gribPrintSec3DP(int *isec0, int *isec3, double *fsec3)
{
  /*

    Print the information in the Bit-Map Section
    (Section 3) of decoded GRIB data.

    Input Parameters:

       isec0  - Array of decoded integers from Section 0

       isec3  - Array of decoded integers from Section 3

       fsec3  - Array of decoded floats from Section 3


    Converted from EMOS routine GRPRS3.

       Uwe Schulzweida   MPIfM   01/04/2001

  */

  UNUSED(isec0);

  grsdef();

  fprintf(grprsm, " \n");
  fprintf(grprsm, " Section 3 - Bit-map Section.\n");
  fprintf(grprsm, " -------------------------------------\n");

  if (isec3[0] != 0)
    fprintf(grprsm, " Predetermined bit-map number.                %9d\n", isec3[0]);
  else
    fprintf(grprsm, " No predetermined bit-map.\n");

  fprintf(grprsm, " Missing data value for integer data.    %14d\n", isec3[1]);

  fprintf(grprsm, " Missing data value for real data. %20.6g\n", fsec3[1]);
}

void
gribPrintSec3SP(int *isec0, int *isec3, float *fsec3sp)
{
  double fsec3[2];

  fsec3[0] = fsec3sp[0];
  fsec3[1] = fsec3sp[1];

  gribPrintSec3DP(isec0, isec3, fsec3);
}

void
gribPrintSec4DP(int *isec0, int *isec4, double *fsec4)
{
  /*

    Print the information in the Binary Data Section
    (Section 4) of decoded GRIB data.

    Input Parameters:

       isec0  - Array of decoded integers from Section 0

       isec4  - Array of decoded integers from Section 4

       fsec4  - Array of decoded floats from Section 4


    Converted from EMOS routine GRPRS4.

       Uwe Schulzweida   MPIfM   01/04/2001

  */
  int inum;

  UNUSED(isec0);

  grsdef();

  /*
    -----------------------------------------------------------------
    Section 1 . Print integer information from isec4.
    -----------------------------------------------------------------
  */
  fprintf(grprsm, " \n");
  fprintf(grprsm, " Section 4 - Binary Data  Section.\n");
  fprintf(grprsm, " -------------------------------------\n");

  fprintf(grprsm, " Number of data values coded/decoded.         %9d\n", isec4[0]);
  fprintf(grprsm, " Number of bits per data value.               %9d\n", isec4[1]);
  fprintf(grprsm, " Type of data       (0=grid pt, 128=spectral).%9d\n", isec4[2]);
  fprintf(grprsm, " Type of packing    (0=simple, 64=complex).   %9d\n", isec4[3]);
  fprintf(grprsm, " Type of data       (0=float, 32=integer).    %9d\n", isec4[4]);
  fprintf(grprsm, " Additional flags   (0=none, 16=present).     %9d\n", isec4[5]);
  fprintf(grprsm, " Reserved.                                    %9d\n", isec4[6]);
  fprintf(grprsm, " Number of values   (0=single, 64=matrix).    %9d\n", isec4[7]);
  fprintf(grprsm, " Secondary bit-maps (0=none, 32=present).     %9d\n", isec4[8]);
  fprintf(grprsm, " Values width       (0=constant, 16=variable).%9d\n", isec4[9]);
  /*
    If complex packing ..
  */
  if (isec4[3] == 64)
    {
      if (isec4[2] == 128)
        {
          fprintf(grprsm, " Byte offset of start of packed data (N).     %9d\n", isec4[15]);
          fprintf(grprsm, " Power (P * 1000).                            %9d\n", isec4[16]);
          fprintf(grprsm, " Pentagonal resolution parameter J for subset.%9d\n", isec4[17]);
          fprintf(grprsm, " Pentagonal resolution parameter K for subset.%9d\n", isec4[18]);
          fprintf(grprsm, " Pentagonal resolution parameter M for subset.%9d\n", isec4[19]);
        }
      else
        {
          fprintf(grprsm, " Bits number of 2nd order values    (none=>0).%9d\n", isec4[10]);
          fprintf(grprsm, " General extend. 2-order packing (0=no,8=yes).%9d\n", isec4[11]);
          fprintf(grprsm, " Boustrophedonic ordering        (0=no,4=yes).%9d\n", isec4[12]);
          fprintf(grprsm, " Spatial differencing order          (0=none).%9d\n", isec4[13] + isec4[14]);
        }
    }
  /*
    Number of non-missing values
  */
  if (isec4[20] != 0) fprintf(grprsm, " Number of non-missing values                 %9d\n", isec4[20]);
  /*
    Information on matrix of values , if present.
  */
  if (isec4[7] == 64)
    {
      fprintf(grprsm, " First dimension (rows) of each matrix.       %9d\n", isec4[49]);
      fprintf(grprsm, " Second dimension (columns) of each matrix.   %9d\n", isec4[50]);
      fprintf(grprsm, " First dimension coordinate values definition.%9d\n", isec4[51]);
      fprintf(grprsm, " (Code Table 12)\n");
      fprintf(grprsm, " NC1 - Number of coefficients for 1st dimension.%7d\n", isec4[52]);
      fprintf(grprsm, " Second dimension coordinate values definition.%8d\n", isec4[53]);
      fprintf(grprsm, " (Code Table 12)\n");
      fprintf(grprsm, " NC2 - Number of coefficients for 2nd dimension.%7d\n", isec4[54]);
      fprintf(grprsm, " 1st dimension physical signifance (Table 13). %8d\n", isec4[55]);
      fprintf(grprsm, " 2nd dimension physical signifance (Table 13).%8d\n", isec4[56]);
    }
  /*
    -----------------------------------------------------------------
    Section 2. Print values from fsec4.
    -----------------------------------------------------------------
  */

  inum = isec4[0];
  if (inum < 0) inum = -inum;
  if (inum > 20) inum = 20;
  /*
    Print first inum values.
  */
  fprintf(grprsm, " \n");
  fprintf(grprsm, " First %4d data values.\n", inum);

  if (isec4[4] == 0)
    {
      /*
        Print real values ...
      */
      for (int j = 0; j < inum; ++j)
        {
          if (fabs(fsec4[j]) > 0)
            {
              if (fabs(fsec4[j]) >= 0.1 && fabs(fsec4[j]) <= 1.e8)
                fprintf(grprsm, " %#16.8G    \n", fsec4[j]);
              else
                fprintf(grprsm, " %#20.8E\n", fsec4[j]);
            }
          else
            fprintf(grprsm, " %#16.0f    \n", fabs(fsec4[j]));
        }
    }
  else
    {
      /*
        Print integer values ...
      */
      fprintf(grprsm, " Print of integer values not supported\n");
      /*
        CALL SETPAR(IBIT,IDUM,IDUM)
        DO 212 J=1,INUM
           INSPT = 0
           CALL INXBIT(IVALUE,1,INSPT,FSEC4(J),1,IBIT,IBIT,'C',IRET)
           WRITE (*,9033) IVALUE
 9033 FORMAT(' ',I15)
  212   CONTINUE
      ENDIF
      */
    }
}

void
gribPrintSec4SP(int *isec0, int *isec4, float *fsec4sp)
{
  double fsec4[20];

  int inum = isec4[0];
  if (inum < 0) inum = -inum;
  if (inum > 20) inum = 20;

  for (int j = 0; j < inum; ++j) fsec4[j] = fsec4sp[j];

  gribPrintSec4DP(isec0, isec4, fsec4);
}

void
gribPrintSec4Wave(int *isec4)
{
  /*

    Print the wave coordinate information in the Binary Data
    Section (Section 4) of decoded GRIB data.

    Input Parameters:

       isec4 - Array of decoded integers from Section 4

    Comments:

       Wave coordinate information held in isec4 are 32-bit floats,
       hence the PTEMP and NTEMP used for printing are 4-byte variables.


    Converted from EMOS routine GRPRS4W.

       Uwe Schulzweida   MPIfM   01/04/2001

  */
  int ntemp[100];
  float *ptemp;

  grsdef();

  /*
    -----------------------------------------------------------------
    Section 1 . Print integer information from isec4.
    -----------------------------------------------------------------
  */
  fprintf(grprsm, " Coefficients defining first dimension coordinates:\n");
  for (int jloop = 0; jloop < isec4[52]; jloop++)
    {
      ntemp[jloop] = isec4[59 + jloop];
      ptemp = (float *) &ntemp[jloop];
      fprintf(grprsm, "%20.10f\n", *ptemp);
    }
  fprintf(grprsm, " Coefficients defining second dimension coordinates:\n");
  for (int jloop = 0; jloop < isec4[54]; jloop++)
    {
      ntemp[jloop] = isec4[59 + isec4[52] + jloop];
      ptemp = (float *) &ntemp[jloop];
      fprintf(grprsm, "%20.10f\n", *ptemp);
    }
}
