external/io_grib1/MEL_grib1/make_grib_log.c

   1 /* Program    : make_grib_log (was printer.c)
   2    Programmer : Todd J. Kienitz, SAIC
   3    Date       : January 10, 1996
   4    Purpose    : To produce the information file output of the GRIB message.
   5    Revisions  :
   6    04/17/96 Steve Lowe, SAIC:  modified data print-out
   7    04/22/96 Alice Nakajima (ATN), SAIC:  added BMS summary
   8    12/12/96 ATN:  implement combined Decoder/Encdoer structs
   9             replaced  (table2 tab2[], table3 tab3[], tables mgotab);
  10    06/14/97 ATN:  print upto encoded pricision.
  11    02/22/98 ATN:  replace projection id with constants, add printing for
  12      prjns: Rotated Lat/Lon, Stretched Lat/Lon, Stretched Rotated Lat/Lon,
  13      Rotated Gaussian, Stretched Gauss, Stretched Rotated Gaussian ,
  14      Oblique Lambert, and for Albers equal-area.
  15    09/10/98 ATN:  extension flag printing.
  16 */
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include "grib_lookup.h"        /* combined encoder/decoder structs */
  21 #include "dprints.h"            /* for debug printing */
  22 #include "gribfuncs.h"          /* prototypes */
  23
  24 /*
  25 **********************************************************************
  26 * A. FUNCTION: make_grib_log
  27 *      Produces debug file GRIB.log from the GRIB message in the Grib Header
  28 *
  29 *    INTERFACE:
  30 *      int   make_grib_log (input_fn, lookup_fn, msg_length, offset,
  31 *                           pds, gds, bds, bms, grib_data, errmsg)
  32 *
  33 *    ARGUMENTS (I=input, O=output, I&O=input and output):
  34 *      (I) char *input_fn;            name of input GRIB file
  35 *      (I) char *lookup_fn;           name of Lookup file, nil if not used
  36 *      (I) unsigned long msg_length;  total length of GRIB message
  37 *      (I) long offset;               starting location of GRIB message in bytes
  38 *      (I) PDS_INPUT pds;             product definition section structure
  39 *      (I) grid_desc_sec gds;         grid description section structure
  40 *      (I) BDS_HEAD_INPUT bds;        binary data section header structure
  41 *      (I) BMS_INPUT bms;             bit map definition section structure
  42 *      (I) float *grib_data;          array of decoded data
  43 *
  44 *     ACCESSES GLOBAL VARS:
  45 *        int UseTables;
  46 *            set to one if lookup table used
  47 *        CTRS_DEFN  db_ctr_tbl[NCTRS];
  48 *            predefined array holding Originating Center info
  49 *        PARM_DEFN db_parm_tbl [MAX_PARM_TBLS * NPARM];
  50 *            predefined arr of Parameter info
  51 *        LVL_DEFN db_lvl_tbl [NLVL];
  52 *            predefined arr of Level info struct
  53 *        MODEL_DEFN db_mdl_tbl [NMODEL];
  54 *            predefined arr of Model info struct
  55 *        GEOM_DEFN db_geom_tbl [NGEOM];
  56 *            predefined arr of Geometry info struct
  57 *
  58 *     RETURN CODE:
  59 *       0> no errors; file GRIB.log has been created;
  60 *       1> error, errmsg filled;
  61 **********************************************************************
  62 */
  63 extern int        UseTables;     /* set means use lookup tbl defns */
  64 extern PARM_DEFN  db_parm_tbl[]; /* parameter conversion info */
  65 extern MODEL_DEFN db_mdl_tbl[];  /* model conversion info */
  66 extern LVL_DEFN   db_lvl_tbl[];  /* level conversion info */
  67 extern GEOM_DEFN  db_geom_tbl[]; /* Geom conversion info  */
  68 extern CTR_DEFN   db_ctr_tbl[];  /* Ctr conversion info  */
  69
  70 #if PROTOTYPE_NEEDED
  71 int   make_grib_log (   char *input_fn,
  72                         char *lookup_fn,
  73                         unsigned long msg_length,
  74                         long offset,
  75                         PDS_INPUT pds,
  76                         grid_desc_sec gds,
  77                         BDS_HEAD_INPUT bds,
  78                         BMS_INPUT bms,
  79                         float *grib_data,
  80                         char *errmsg)
  81 #else
  82 int   make_grib_log (input_fn, lookup_fn, msg_length, offset,
  83                      pds, gds, bds, bms, grib_data,errmsg)
  84         char *input_fn;
  85         char *lookup_fn;
  86         unsigned long msg_length;
  87         long offset;
  88         PDS_INPUT pds;
  89         grid_desc_sec gds;
  90         BDS_HEAD_INPUT bds;
  91         BMS_INPUT bms;
  92         float *grib_data;
  93         char *errmsg;
  94 #endif
  95
  96 {
  97   char *func="make_grib_log";
  98   int i, indx, k, fd, numpts=100;
  99   float dsf, res, min, max;
 100   FILE *fp;
 101
 102 /*
 103 *
 104 * A.0   DEBUG printing
 105 */
 106   DPRINT1 ("Entering %s\n", func);
 107
 108 /*
 109 *
 110 * A.1   OPEN file "GRIB.log" in APPEND mode
 111 */
 112   fp=fopen ("GRIB.log", "a+");
 113   if (!fp) {
 114         DPRINT1("%s:  failed to open 'GRIB.log' for appending, skip logfile\n",
 115         func);
 116         sprintf (errmsg, "%s:  failed to open 'GRIB.log'\n", func);
 117         return (1);
 118         }
 119
 120 /*
 121 *
 122 * A.2   WRITE Indicator Section information to file
 123 *       !message length
 124 *       !GRIB Edition number
 125 */
 126   fseek(fp, 0L, 2);
 127   if (ftell(fp) == 0L)
 128   fprintf (fp,  "%s: InFile= %s\n%s: Lookup=%d, fn='%s'\n\n" ,
 129            func, input_fn, func,UseTables, lookup_fn);
 130
 131   fprintf (fp, "**** VALID MESSAGE FOUND AT %ld BYTES ****\n" , offset);
 132
 133   fprintf(fp, "\n********* SECTION 0 IDS *********\n" );
 134   fprintf(fp, "Total Message length = %ld\nEdition Number = %d\n",
 135            msg_length, pds.usEd_num);
 136 /*
 137 *
 138 * A.3   WRITE Product Definition Section information to file
 139 *       !Section length
 140 *       !Parameter Table version
 141 *       !Parameter Sub-Table version if defined and flagged by Extension flag
 142 *       !Tracking id if defined and flagged by Extension flag
 143 */
 144   fprintf(fp, "\n********* SECTION 1 PDS *********\n" \
 145         "Section length = %d\nTable version = %d\n",
 146          pds.uslength, pds.usParm_tbl);
 147
 148   if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG )
 149    {
 150     if (pds.usSub_tbl != 0)
 151         fprintf(fp,"Local Table version = %d\n",pds.usSub_tbl);
 152     if(pds.usTrack_num  != 0)
 153         fprintf(fp,"Tracking ID = %d\n",pds.usTrack_num);
 154    }
 155
 156 /*
 157 *       !Originating Center id
 158 *       !IF (using tables) Name of Originating Center
 159 */
 160   fprintf(fp,"Originating Center id = %d\n",pds.usCenter_id);
 161   if (UseTables)
 162      if ( db_ctr_tbl[pds.usCenter_id].ctr_dsc[0] )
 163         fprintf(fp,"Originating Center = %s\n",
 164         db_ctr_tbl[pds.usCenter_id].ctr_dsc);
 165      else
 166         fprintf(fp,"Originating Center ID %d not defined in current table.\n",
 167         pds.usCenter_id);
 168
 169 /*
 170 *       !Sub-Table Entry for Originating Center if non-zero and if
 171 *       !extension flag is set
 172 */
 173   if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG &&
 174         pds.usCenter_sub != 0)
 175      fprintf(fp,"Sub-Table Entry Originating Center = %d\n",pds.usCenter_sub);
 176
 177 /*
 178 *       !Extension flag
 179 */
 180   fprintf(fp,"Extension flag = %d (extensions %s)\n", pds.usExt_flag,
 181   (pds.usExt_flag == (unsigned short)EXTENSION_FLAG ? "used" : "not used"));
 182
 183 /*
 184 *       !Model Identification
 185 *       !IF (using tables) Model Description
 186 */
 187   fprintf(fp,"Model id = %d\n",pds.usProc_id);
 188   if (UseTables)
 189      if ( db_mdl_tbl[pds.usProc_id].grib_dsc[0] )
 190         fprintf(fp,"Model Description = %s\n",
 191         db_mdl_tbl[pds.usProc_id].grib_dsc);
 192      else
 193         fprintf(fp,"Model ID %d not defined in current table.\n",
 194         pds.usProc_id);
 195
 196 /*
 197 *       !Grid Identification
 198 *       !IF (using tables) Grid Description
 199 */
 200   fprintf(fp,"Grid id = %d\n",pds.usGrid_id);
 201   if (UseTables)
 202      if ( db_geom_tbl[pds.usGrid_id].grib_dsc[0] )
 203         fprintf(fp,"Grid Description = %s\n",
 204         db_geom_tbl[pds.usGrid_id].grib_dsc);
 205      else
 206         fprintf(fp,"Grid ID %d not defined in current table.\n",
 207         pds.usGrid_id);
 208
 209 /*
 210 *       !Parameter Identification
 211 */
 212   fprintf(fp,"Parameter id = %d\n",pds.usParm_id);
 213
 214 /*
 215 *       !IF (usExt_flag is set AND
 216 *       !    (Parm id between 250 and 254) AND (Sub Parm ID defined)))
 217 *       !    PRINT Parm_sub
 218 *       !ENDIF
 219 */
 220   if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG &&
 221           pds.usParm_id>=250 && pds.usParm_id<=254 && pds.usParm_sub!=0)
 222           fprintf(fp,"Parameter sub-id = %d\n",pds.usParm_sub);
 223
 224
 225 /*
 226 *       !IF (using lookup table) THEN
 227 *       !  CALCULATE index in Parm Conversion Array to use
 228 *       !  let index= (usParm_Id - 249)*256 + usParm_sub;
 229 *       !
 230 *       !   IF this index in Parm Conversion Array is defined THEN
 231 *       !       PRINT its grib_dsc and grib_unit_dsc
 232 *       !   ELSE
 233 *       !       PRINT it's not defined mesage
 234 *       !   ENDIF
 235 *       !ENDIF
 236 */
 237   if(UseTables)
 238    {
 239       indx = PARMTBL_INDX (pds.usParm_id, pds.usParm_sub);
 240
 241       if ( db_parm_tbl[indx].grib_dsc[0] ) {
 242          fprintf(fp,"Parameter name = %s\n",db_parm_tbl[indx].grib_dsc);
 243          fprintf(fp,"Parameter units = %s\n",db_parm_tbl[indx].grib_unit_dsc);
 244           }
 245       else fprintf(fp,"Parameter ID %d not defined in current table.\n",
 246          pds.usParm_id);
 247     }
 248
 249 /*
 250 *       !Level Id
 251 *       !IF (using tables)
 252 *       !  Level description
 253 *       !  SWITCH (number of octets to store Height1)
 254 *       !     2: Level = Height1
 255 *       !     1: Bottom of Layer = Height1
 256 *       !        Top of Layer = Height2
 257 *       !     0: (no Height value required)
 258 *       !     default: (corrupt table entry or message)
 259 *       !  ENDSWITCH
 260 *       !ELSE (not using tables)
 261 *       !  Level = Height1  (Level assumed)
 262 *       !ENDIF
 263 */
 264   fprintf(fp,"Level_type = %d\n",pds.usLevel_id);
 265   if(UseTables) {
 266     if ( db_lvl_tbl[pds.usLevel_id].grib_dsc[0] ) {
 267        fprintf(fp,"Level description = %s\n",
 268                 db_lvl_tbl[pds.usLevel_id].grib_dsc);
 269        switch(db_lvl_tbl[pds.usLevel_id].num_octets){
 270          case 2:
 271            fprintf(fp,"%s = %u\n",
 272            db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1);
 273            break;
 274          case 1:
 275            fprintf(fp,"%s = %u\n%s = %u\n",
 276            db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1,
 277            db_lvl_tbl[pds.usLevel_id].lvl_name_2, pds.usHeight2);
 278            break;
 279          case 0:
 280            break;
 281          default:
 282            fprintf(fp,"***Number of octets for table 3 undefined - possibly "
 283                    "corrupt dataset.***\n");
 284        }
 285     }else
 286        fprintf(fp,"Level ID %d not defined in current table.\n",
 287         pds.usLevel_id);
 288   } /* end UseTables 'if' statement */
 289   else fprintf(fp,"Level = %u\n",pds.usHeight1);
 290
 291 /*
 292 *       !Reference Date/Time:
 293 *       !  Century
 294 *       !  Year
 295 *       !  Month
 296 *       !  Day
 297 *       !  Hour
 298 *       !  Minute
 299 *       !  Second if defined
 300 */
 301   fprintf(fp,
 302      "Reference Date/Time of Data Set:\n" \
 303      "   Century = %d\n   Year = %d\n   Month = %d\n   Day = %d\n"\
 304      "   Hour = %d\n   Minute = %d\n",
 305      pds.usCentury,pds.usYear,pds.usMonth,pds.usDay,pds.usHour,pds.usMinute);
 306
 307   if(pds.usExt_flag == (unsigned short)EXTENSION_FLAG)
 308         fprintf(fp,"   Second = %d\n",pds.usSecond);
 309
 310 /*
 311 *       !Forecast Time Unit
 312 *       !  Forecast Period 1
 313 *       !  Forecast Period 2
 314 */
 315   switch(pds.usFcst_unit_id){
 316     case 0: fprintf(fp,"Forecast Time Unit = Minute\n"); break;
 317     case 1: fprintf(fp,"Forecast Time Unit = Hour\n"); break;
 318     case 2: fprintf(fp,"Forecast Time Unit = Day\n"); break;
 319     case 3: fprintf(fp,"Forecast Time Unit = Month\n"); break;
 320     case 4: fprintf(fp,"Forecast Time Unit = Year\n"); break;
 321     case 5: fprintf(fp,"Forecast Time Unit = Decade (10 years)\n"); break;
 322     case 6: fprintf(fp,"Forecast Time Unit = Normal (30 years)\n"); break;
 323     case 7: fprintf(fp,"Forecast Time Unit = Century (100 years)\n"); break;
 324     case 254: fprintf(fp,"Forecast Time Unit = Second\n"); break;
 325     default: fprintf(fp,"Forecast Time Unit = UNDEFINED!!\n");
 326   }
 327   fprintf(fp,"   Forecast Period 1 = %d\n",pds.usP1);
 328   fprintf(fp,"   Forecast Period 2 = %d\n",pds.usP2);
 329
 330 /*
 331 *       !Time Range Indicator
 332 *       !Number in Average
 333 *       !Number Missing
 334 */
 335   fprintf(fp,"Time Range = %d\n",pds.usTime_range);
 336   fprintf(fp,"Number in Average = %d\n",pds.usTime_range_avg);
 337   fprintf(fp,"Number Missing = %d\n",pds.usTime_range_mis);
 338
 339 /*
 340 *       !Decimal Scale Factor
 341 */
 342   fprintf(fp,"Decimal Scale Factor = %d\n",pds.sDec_sc_fctr);
 343
 344 /*
 345 *
 346 * A.4   IF (GDS included) THEN
 347 * A.4.1    WRITE Grid Definition Section information to file
 348 *            !Section length
 349 *            !Parm_nv
 350 *            !Parm_pv_pl
 351 *            !Data type
 352 */
 353   if(pds.usGds_bms_id >> 7 & 1) {
 354
 355      fprintf(fp,"\n********* SECTION 2 GDS *********\n");
 356      fprintf(fp,"Section length = %d\n",gds.head.uslength);
 357      fprintf(fp,"Parm_nv = %d\n",gds.head.usNum_v);
 358      fprintf(fp,"Parm_pv_pl = %d\n",gds.head.usPl_Pv);
 359      fprintf(fp,"Data_type = %d\n",gds.head.usData_type);
 360
 361 /*
 362 * A.4.2    SWITCH (Data Type, Table 6)
 363 *                 !  For each Data Type, write the following to file:
 364 *                 !     Number of points along rows/columns of grid
 365 *                 !     Reference Lat/Lon information
 366 *                 !     Resolution and Component Flags (Table 7)
 367 *                 !       Direction increments if given
 368 *                 !       Assumption of Earth shape
 369 *                 !       U&V component orientation
 370 *                 !     Scanning mode flags (Table 8)
 371 *              Default: Projection not supported, exit;
 372 */
 373      fprintf(fp,"Projection = %s\n", prjn_name[gds.head.usData_type]);
 374      switch(gds.head.usData_type)
 375      {
 376
 377 /*
 378 *               Case  0: Lat/Lon projection
 379 *               Case 10: Rotated Lat/Lon projection
 380 *               Case 20: Stretched Lat/Lon projection
 381 *               Case 30: Stretched Rotated Lat/Lon projection
 382 */
 383     case LATLON_PRJ:            /* Lat/Lon Grid */
 384     case ROT_LATLON_PRJ:        /* Rotated Lat/Lon */
 385     case STR_LATLON_PRJ:        /* Stretched Lat/Lon */
 386     case STR_ROT_LATLON_PRJ :  /* Stretched and Rotated Lat/Lon */
 387
 388       fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi);
 389       fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj);
 390       fprintf(fp,"Latitude of first grid point = %.3f deg\n",
 391       ((float)gds.llg.lLat1)/1000.);
 392       fprintf(fp,"Longitude of first grid point = %.3f deg\n",
 393       ((float)gds.llg.lLon1)/1000.);
 394       fprintf(fp,"Latitude of last grid point = %.3f deg\n",
 395       ((float)gds.llg.lLat2)/1000.);
 396       fprintf(fp,"Longitude of last grid point = %.3f deg\n",
 397       ((float)gds.llg.lLon2)/1000.);
 398
 399       fprintf(fp,"Resolution and Component Flags: \n");
 400       if ((gds.llg.usRes_flag >> 7) & 1) {
 401          fprintf(fp,"    Longitudinal increment = %f deg\n",
 402         ((float)gds.llg.iDi)/1000.);
 403          fprintf(fp,"    Latitudinal increment = %f deg\n",
 404         ((float)gds.llg.iDj)/1000.);
 405       }else fprintf(fp,"    Direction increments not given.\n");
 406       if ((gds.llg.usRes_flag >> 6) & 1)
 407            fprintf(fp,"    Earth assumed oblate spherical.\n");
 408       else fprintf(fp,"    Earth assumed spherical.\n");
 409       if ((gds.llg.usRes_flag >> 3) & 1)
 410            fprintf(fp,"    U&V components resolved relative to +I and "
 411                    "+J\n");
 412       else fprintf(fp,"    U&V components resolved relative to east "
 413                    "and north.\n");
 414
 415       fprintf(fp,"Scanning Mode Flags: \n");
 416       if ((gds.llg.usScan_mode >> 7) & 1)
 417            fprintf(fp,"    Points scan in -I direction.\n");
 418       else fprintf(fp,"    Points scan in +I direction.\n");
 419       if ((gds.llg.usScan_mode >> 6) & 1)
 420            fprintf(fp,"    Points scan in +J direction.\n");
 421       else fprintf(fp,"    Points scan in -J direction.\n");
 422       if ((gds.llg.usScan_mode >> 5) & 1)
 423            fprintf(fp,"    Adjacent points in J direction are "
 424                    "consecutive.\n");
 425       else fprintf(fp,"    Adjacent points in I direction are "
 426                    "consecutive.\n");
 427
 428       /* added 02/98
 429         This code pertains only to the Stretch/Rotate grids,
 430         so skip over if it's a LATLON_PRJN type;
 431         */
 432       if (gds.head.usData_type != LATLON_PRJ)
 433         {
 434         fprintf(fp,"    Latitude of southern pole = %.3f deg\n",
 435         ((float)gds.llg.lLat_southpole)/1000.);
 436         fprintf(fp,"    Longitude of southern pole = %.3f deg\n",
 437         ((float)gds.llg.lLon_southpole)/1000.);
 438
 439         /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
 440            a single precision floating point value */
 441         fprintf(fp,"    Angle of rotation = %.3f\n",
 442         grib_ibm_local ((unsigned long) gds.llg.lRotate));
 443
 444         fprintf(fp,"    Latitude of pole of stretching = %.3f deg\n",
 445         ((float)gds.llg.lPole_lat)/1000.);
 446         fprintf(fp,"    Longitude of pole of stretching = %.3f deg\n",
 447         ((float)gds.llg.lPole_lon)/1000.);
 448
 449         /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
 450            a single precision floating point value */
 451         fprintf(fp,"    Stretching factor = %.3f\n",
 452         grib_ibm_local ((unsigned long)gds.llg.lStretch));
 453         }
 454       break;
 455
 456 /*
 457 *               Case 1: Mercator Projection
 458 */
 459     case MERC_PRJ:    /* Mercator Projection Grid    */
 460
 461       fprintf(fp,"Number of points along a parallel = %d\n",gds.merc.cols);
 462       fprintf(fp,"Number of points along a meridian = %d\n",gds.merc.rows);
 463       fprintf(fp,"Latitude of first grid point = %.3f deg\n",
 464       ((float)gds.merc.first_lat)/1000.);
 465       fprintf(fp,"Longitude of first grid point = %.3f deg\n",
 466       ((float)gds.merc.first_lon)/1000.);
 467       fprintf(fp,"Latitude of last grid point = %.3f deg\n",
 468       ((float)gds.merc.La2)/1000.);
 469       fprintf(fp,"Longitude of last grid point = %.3f deg\n",
 470       ((float)gds.merc.Lo2)/1000.);
 471       fprintf(fp,"Latitude of intersection with Earth = %.3f deg\n",
 472       ((float)gds.merc.latin)/1000.);
 473
 474       fprintf(fp,"Resolution and Component Flags: \n");
 475       if ((gds.merc.usRes_flag >> 7) & 1) {
 476          fprintf(fp,"    Longitudinal increment = %f deg\n",
 477          ((float)gds.merc.lon_inc)/1000.);
 478          fprintf(fp,"    Latitudinal increment = %f deg\n",
 479          ((float)gds.merc.lat_inc)/1000.);
 480       }else fprintf(fp,"    Direction increments not given.\n");
 481       if ((gds.merc.usRes_flag >> 6) & 1)
 482            fprintf(fp,"    Earth assumed oblate spherical.\n");
 483       else fprintf(fp,"    Earth assumed spherical.\n");
 484       if ((gds.merc.usRes_flag >> 3) & 1)
 485            fprintf(fp,"    U&V components resolved relative to +I and "
 486                    "+J\n");
 487       else fprintf(fp,"    U&V components resolved relative to east "
 488                    "and north.\n");
 489
 490       fprintf(fp,"Scanning Mode Flags: \n");
 491       if ((gds.merc.usScan_mode >> 7) & 1)
 492            fprintf(fp,"    Points scan in -I direction.\n");
 493       else fprintf(fp,"    Points scan in +I direction.\n");
 494       if ((gds.merc.usScan_mode >> 6) & 1)
 495            fprintf(fp,"    Points scan in +J direction.\n");
 496       else fprintf(fp,"    Points scan in -J direction.\n");
 497       if ((gds.merc.usScan_mode >> 5) & 1)
 498            fprintf(fp,"    Adjacent points in J direction are "
 499                    "consecutive.\n");
 500       else fprintf(fp,"    Adjacent points in I direction are "
 501                    "consecutive.\n");
 502       break;
 503
 504 /*
 505 *               Case  3: Lambert Conformal Projection
 506 *               Case 13: Oblique Lambert Conformal Projection
 507 *               Case  8: Alberts equal-area secant/tangent conic/bipolar Prj
 508 */
 509     case LAMB_PRJ:              /* Lambert Conformal */
 510     case OBLIQ_LAMB_PRJ:        /* Oblique Lambert Conformal */
 511     case ALBERS_PRJ:            /* Albers equal-area */
 512
 513       fprintf(fp,"Number of points along X-axis = %d\n",gds.lam.iNx);
 514       fprintf(fp,"Number of points along Y-axis = %d\n",gds.lam.iNy);
 515       fprintf(fp,"Latitude of first grid point = %.3f deg\n",
 516       ((float)gds.lam.lLat1)/1000.);
 517       fprintf(fp,"Longitude of first grid point = %.3f deg\n",
 518       ((float)gds.lam.lLon1)/1000.);
 519       fprintf(fp,"Orientation of grid = %.3f deg\n",
 520       ((float)gds.lam.lLon_orient)/1000.);
 521       fprintf(fp,"First Latitude Cut = %.3f deg\n",
 522       ((float)gds.lam.lLat_cut1)/1000.);
 523       fprintf(fp,"Second Latitude Cut = %.3f deg\n",
 524       ((float)gds.lam.lLat_cut2)/1000.);
 525
 526       fprintf(fp,"Resolution and Component Flags: \n");
 527       if ((gds.lam.usRes_flag >> 7) & 1) {
 528             fprintf(fp,"    X-direction increment = %d meters\n",
 529             gds.lam.ulDx);
 530             fprintf(fp,"    Y-direction increment = %d meters\n",
 531             gds.lam.ulDy);
 532       }else fprintf(fp,"    Direction increments not given.\n");
 533       if ((gds.lam.usRes_flag >> 6) & 1)
 534            fprintf(fp,"    Earth assumed oblate spherical.\n");
 535       else fprintf(fp,"    Earth assumed spherical.\n");
 536       if ((gds.lam.usRes_flag >> 3) & 1)
 537            fprintf(fp,"    U&V components resolved relative to +I and "
 538                    "+J\n");
 539       else fprintf(fp,"    U&V components resolved relative to east "
 540                    "and north.\n");
 541
 542       fprintf(fp,"Scanning Mode Flags: \n");
 543       if ((gds.lam.usScan_mode >> 7) & 1)
 544            fprintf(fp,"    Points scan in -I direction.\n");
 545       else fprintf(fp,"    Points scan in +I direction.\n");
 546       if ((gds.lam.usScan_mode >> 6) & 1)
 547            fprintf(fp,"    Points scan in +J direction.\n");
 548       else fprintf(fp,"    Points scan in -J direction.\n");
 549       if ((gds.lam.usScan_mode >> 5) & 1)
 550            fprintf(fp,"    Adjacent points in J direction are "
 551                    "consecutive.\n");
 552       else fprintf(fp,"    Adjacent points in I direction are "
 553                    "consecutive.\n");
 554
 555       /* 02/98 This code pertains only to the Albers projection */
 556       if (gds.head.usData_type == ALBERS_PRJ)
 557         {
 558         fprintf(fp,"    Latitude of the southern pole = %.3f\n",
 559         ((float)gds.lam.lLat_southpole)/1000.);
 560         fprintf(fp,"    Longitude of the southern pole = %.3f\n",
 561         ((float)gds.lam.lLon_southpole)/1000.);
 562         }
 563       break;
 564
 565 /*
 566 *               Case  4: Gaussian Lat/Lon Projection
 567 *               Case 14: Rotated Gaussian Lat/Lon Projection
 568 *               Case 24: Stretched Gaussian Lat/Lon Projection
 569 *               Case 34: Stretched Rotated Gaussian Lat/Lon Projection
 570 */
 571     case GAUSS_PRJ:             /* Gaussian Latitude/Longitude Grid */
 572     case ROT_GAUSS_PRJ:         /* Rotated Gaussian */
 573     case STR_GAUSS_PRJ :        /* Stretched Gaussian */
 574     case STR_ROT_GAUSS_PRJ :    /* Stretched and Rotated Gaussian */
 575
 576       fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi);
 577       fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj);
 578       fprintf(fp,"Latitude of first grid point = %.3f deg\n",
 579       ((float)gds.llg.lLat1)/1000.);
 580       fprintf(fp,"Longitude of first grid point = %.3f deg\n",
 581       ((float)gds.llg.lLon1)/1000.);
 582       fprintf(fp,"Latitude of last grid point = %.3f deg\n",
 583       ((float)gds.llg.lLat2)/1000.);
 584       fprintf(fp,"Longitude of last grid point = %.3f deg\n",
 585       ((float)gds.llg.lLon2)/1000.);
 586
 587       fprintf(fp,"Resolution and Component Flags: \n");
 588       if ((gds.llg.usRes_flag >> 7) & 1) {
 589          fprintf(fp,"    i direction increment = %f deg\n",
 590         ((float)gds.llg.iDi)/1000.);
 591          fprintf(fp,
 592         "    Number of parallels between pole and equator = %d\n",
 593         gds.llg.iDj);
 594       }else fprintf(fp,"    Direction increments not given.\n");
 595       if ((gds.llg.usRes_flag >> 6) & 1)
 596            fprintf(fp,"    Earth assumed oblate spherical.\n");
 597       else fprintf(fp,"    Earth assumed spherical.\n");
 598       if ((gds.llg.usRes_flag >> 3) & 1)
 599            fprintf(fp,"    U&V components resolved relative to +I and "
 600                    "+J\n");
 601       else fprintf(fp,"    U&V components resolved relative to east "
 602                    "and north.\n");
 603
 604       fprintf(fp,"Scanning Mode Flags: \n");
 605       if ((gds.llg.usScan_mode >> 7) & 1)
 606            fprintf(fp,"    Points scan in -I direction.\n");
 607       else fprintf(fp,"    Points scan in +I direction.\n");
 608       if ((gds.llg.usScan_mode >> 6) & 1)
 609            fprintf(fp,"    Points scan in +J direction.\n");
 610       else fprintf(fp,"    Points scan in -J direction.\n");
 611       if ((gds.llg.usScan_mode >> 5) & 1)
 612            fprintf(fp,"    Adjacent points in J direction are "
 613                    "consecutive.\n");
 614       else fprintf(fp,"    Adjacent points in I direction are "
 615                    "consecutive.\n");
 616
 617       /* added 02/98
 618         This code pertains only to the Stretch/Rotate grids
 619         */
 620       if (gds.head.usData_type != GAUSS_PRJ)
 621         {
 622         fprintf(fp,"    Latitude of southern pole = %.3f deg\n",
 623         ((float)gds.llg.lLat_southpole)/1000.);
 624         fprintf(fp,"    Longitude of southern pole = %.3f deg\n",
 625         ((float)gds.llg.lLon_southpole)/1000.);
 626
 627        /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
 628            a single precision floating point value */
 629         fprintf(fp,"    Angle of rotation = %.3f\n",
 630         grib_ibm_local ((unsigned long) gds.llg.lRotate));
 631
 632         fprintf(fp,"    Latitude of pole of stretching = %.3f deg\n",
 633         (float)gds.llg.lPole_lat);
 634         fprintf(fp,"    Longitude of pole of stretching = %.3f deg\n",
 635         (float)gds.llg.lPole_lon);
 636
 637         /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
 638            a single precision floating point value */
 639         fprintf(fp,"    Stretching factor = %.3f\n",
 640         grib_ibm_local ((unsigned long)gds.llg.lStretch));
 641         }
 642       break;
 643
 644 /*
 645 *               Case 5: Polar Sterographic Projection
 646 */
 647     case POLAR_PRJ:    /* Polar Stereographic Projection Grid    */
 648       fprintf(fp,"Number of points along X-axis = %d\n",gds.pol.usNx);
 649       fprintf(fp,"Number of points along Y-axis = %d\n",gds.pol.usNy);
 650       fprintf(fp,"Latitude of first grid point = %.3f deg\n",
 651         ((float)gds.pol.lLat1)/1000.);
 652       fprintf(fp,"Longitude of first grid point = %.3f deg\n",
 653         ((float)gds.pol.lLon1)/1000.);
 654       fprintf(fp,"Orientation of grid = %.3f deg\n",
 655         ((float)gds.pol.lLon_orient)/1000.);
 656       fprintf(fp,"Projection Center: ");
 657       if ((gds.pol.usProj_flag >> 7) & 1)
 658            fprintf(fp,"South Pole\n");
 659       else fprintf(fp,"North Pole\n");
 660
 661       fprintf(fp,"Resolution and Component Flags: \n");
 662       if ((gds.pol.usRes_flag >> 7) & 1) {
 663          fprintf(fp,"    X-direction grid length = %d meters\n",gds.pol.ulDx);
 664          fprintf(fp,"    Y-direction grid length = %d meters\n",gds.pol.ulDy);
 665       }else fprintf(fp,"    Direction increments not given.\n");
 666       if ((gds.pol.usRes_flag >> 6) & 1)
 667            fprintf(fp,"    Earth assumed oblate spherical.\n");
 668       else fprintf(fp,"    Earth assumed spherical.\n");
 669       if ((gds.pol.usRes_flag >> 3) & 1)
 670            fprintf(fp,"    U&V components resolved relative to +I and "
 671                    "+J\n");
 672       else fprintf(fp,"    U&V components resolved relative to east "
 673                    "and north.\n");
 674
 675       fprintf(fp,"Scanning Mode Flags: \n");
 676       if ((gds.pol.usScan_mode >> 7) & 1)
 677            fprintf(fp,"    Points scan in -I direction.\n");
 678       else fprintf(fp,"    Points scan in +I direction.\n");
 679       if ((gds.pol.usScan_mode >> 6) & 1)
 680            fprintf(fp,"    Points scan in +J direction.\n");
 681       else fprintf(fp,"    Points scan in -J direction.\n");
 682       if ((gds.pol.usScan_mode >> 5) & 1)
 683            fprintf(fp,"    Adjacent points in J direction are "
 684                    "consecutive.\n");
 685       else fprintf(fp,"    Adjacent points in I direction are "
 686                    "consecutive.\n");
 687       break;
 688
 689     default:   /* Bad projection:  ignore & continue */
 690       fprintf(stdout, "\n\n***%s WARNING ***:\nProjection %d is INVALID;\n\n",
 691       func, gds.head.usData_type);
 692
 693       fprintf(fp,"================================================\n"\
 694                  "%s:  projection %d is not currently implemented\n"\
 695                 "================================================\n",
 696       func, gds.head.usData_type);
 697       break;
 698
 699 /*
 700 * A.4.2    ENDSWITCH (Data Type)
 701 */
 702      } /* Switch */
 703
 704   } /* gds included */
 705 /*
 706 *
 707 * A.4     ELSE
 708 *             PRINT no Gds message
 709 * A.4     ENDIF
 710 */
 711   else fprintf(fp,"\n******* NO SECTION 2 GDS *********\n" );
 712
 713
 714 /*
 715 *
 716 * A.5   IF (Bitmap Section is present)
 717 *       THEN
 718 *          WRITE Bitmap Section information to file
 719 *       ELSE
 720 *          PRINT no bms mesg
 721 *       ENDIF
 722 */
 723   if(pds.usGds_bms_id >> 6 & 1) {
 724     fprintf(fp,"\n********* SECTION 3 BMS **********\n" );
 725     fprintf(fp,"Section length = %ld\n", bms.uslength);
 726     if (bms.uslength <= 6)
 727       fprintf(fp,"Bitmap is predefined (Not in message).\n");
 728     else fprintf(fp,"Bitmap is included with message.\n");
 729     fprintf(fp,"Bitmap ID = %d \n", bms.usBMS_id);
 730     fprintf(fp,"Number of unused bits = %d\n", bms.usUnused_bits);
 731     fprintf(fp,"Number of datapoints set = %ld\n", bms.ulbits_set);
 732   }else{
 733     fprintf(fp,"\n******* NO SECTION 3 BMS *********\n" );
 734   }
 735
 736 /*
 737 *
 738 * A.6   WRITE out Binary Data Section Information to file
 739 *       !Section Length
 740 */
 741   fprintf(fp,"\n********* SECTION 4 BDS *********\n" );
 742   fprintf(fp,"Section length = %ld\n",bds.length);
 743
 744 /*
 745 *       !Table 11 Flags
 746 */
 747   fprintf(fp,"Table 11 Flags:\n");
 748   if ((bds.usBDS_flag >> 7) & 1)
 749        fprintf(fp,"    Spherical harmonic coefficients.\n");
 750   else fprintf(fp,"    Grid-point data.\n");
 751   if ((bds.usBDS_flag >> 6) & 1)
 752        fprintf(fp,"    Second-order packing.\n");
 753   else fprintf(fp,"    Simple Packing.\n");
 754   if ((bds.usBDS_flag >> 5) & 1)
 755        fprintf(fp,"    Integer values.\n");
 756   else fprintf(fp,"    Floating point values.\n");
 757   if ((bds.usBDS_flag >> 4) & 1)
 758        fprintf(fp,"    Octet 14 contains additional flag bits.\n");
 759   else fprintf(fp,"    No additional flags at octet 14.\n");
 760
 761 /*
 762 *       !Decimal Scale Factor (Repeated from PDS)
 763 */
 764   fprintf(fp,"\nDecimal Scale Factor = %d\n",pds.sDec_sc_fctr);
 765
 766 /*
 767 *       !Binary Scale Factor
 768 *       !Bit Width
 769 *       !Number of Data Points
 770 */
 771   fprintf(fp,"Binary scale factor = %d\n", bds.Bin_sc_fctr);
 772   fprintf(fp,"Bit width = %d\n", bds.usBit_pack_num);
 773   fprintf(fp,"Number of data points = %ld\n",bds.ulGrid_size);
 774
 775 /*
 776 * A.6.1   WRITE Data Summary to file
 777 *         !Compute Data Min/Max and Resolution
 778 */
 779   dsf = (float) pow( (double) 10, (double) pds.sDec_sc_fctr);
 780   res = (float) pow((double)2,(double)bds.Bin_sc_fctr) / dsf;
 781   min = bds.fReference / dsf;
 782   max = (float) (pow((double)2, (double)bds.usBit_pack_num) - 1);
 783   max = min + max * res;
 784   fprintf(fp,"Data Minimum = %f\n", min );
 785   fprintf(fp,"Data Maximum = %f\n", max );
 786   fprintf(fp,"Resolution = %f\n",res );
 787
 788 /*
 789 *         !Compute Format Specifier for printing Data
 790 */
 791   fd = (int) -1 * (float) log10((double) res) + .5;
 792   if (fd <= 0)
 793   {
 794     fd = 0;
 795     fprintf(fp,"DATA will be displayed as integers (res > 0.1).\n");
 796   }
 797
 798 /*
 799 *         !WRITE First 100 Data Points to file up to Encoded Precision
 800 */
 801   if (bds.ulGrid_size > 1) {
 802   if (bds.ulGrid_size < 100) numpts = bds.ulGrid_size;
 803   fprintf(fp,"\nDATA ARRAY: (first %d)\n",numpts);
 804   if (fd > 0)  {
 805         for (i=0; i<numpts; i=i+5)
 806           if (i + 5 < numpts)
 807               fprintf(fp, "%03d-  %.*f  %.*f  %.*f  %.*f  %.*f\n",
 808               i,fd,grib_data[i],fd,grib_data[i+1],fd,
 809               grib_data[i+2],fd,grib_data[i+3],fd,grib_data[i+4] );
 810           else  {
 811               fprintf(fp,"%03d-", i);
 812               while (i < numpts) fprintf(fp,"  %.*f", fd, grib_data[i++]);
 813               fprintf(fp,"\n");
 814             }
 815         }
 816   else  {
 817         for (i=0; i<numpts; i=i+5)
 818           if (i + 5 < numpts)
 819               fprintf(fp, "%03d-  %.0f  %.0f  %.0f  %.0f  %.0f\n",
 820               i,grib_data[i],grib_data[i+1],
 821               grib_data[i+2],grib_data[i+3],grib_data[i+4] );
 822           else  {
 823               fprintf(fp,"%03d-", i);
 824               while (i < numpts) fprintf(fp,"  %.0f", grib_data[i++]);
 825               fprintf(fp,"\n");
 826             }
 827         }
 828
 829   }
 830
 831   fprintf (fp,"\n******** END OF MESSAGE *******\n\n");
 832
 833 /*
 834 *
 835 * A.7   CLOSE file
 836 */
 837   fclose (fp);
 838
 839 /*
 840 *
 841 * A.8   DEBUG printing
 842 */
 843   DPRINT1 ("Leaving  %s, no errors\n", func);
 844   return (0);
 845
 846 /*
 847 * END OF FUNCTION
 848 *
 849 *
 850 */
 851 }