/wrfv2_fire/external/io_grib1/MEL_grib1/make_grib_log.c
C | 851 lines | 505 code | 74 blank | 272 comment | 121 complexity | 7073bb2d627ee26c16e5511daf088e6f MD5 | raw file
Possible License(s): AGPL-1.0
- /* Program : make_grib_log (was printer.c)
- Programmer : Todd J. Kienitz, SAIC
- Date : January 10, 1996
- Purpose : To produce the information file output of the GRIB message.
- Revisions :
- 04/17/96 Steve Lowe, SAIC: modified data print-out
- 04/22/96 Alice Nakajima (ATN), SAIC: added BMS summary
- 12/12/96 ATN: implement combined Decoder/Encdoer structs
- replaced (table2 tab2[], table3 tab3[], tables mgotab);
- 06/14/97 ATN: print upto encoded pricision.
- 02/22/98 ATN: replace projection id with constants, add printing for
- prjns: Rotated Lat/Lon, Stretched Lat/Lon, Stretched Rotated Lat/Lon,
- Rotated Gaussian, Stretched Gauss, Stretched Rotated Gaussian ,
- Oblique Lambert, and for Albers equal-area.
- 09/10/98 ATN: extension flag printing.
- */
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #include "grib_lookup.h" /* combined encoder/decoder structs */
- #include "dprints.h" /* for debug printing */
- #include "gribfuncs.h" /* prototypes */
- /*
- **********************************************************************
- * A. FUNCTION: make_grib_log
- * Produces debug file GRIB.log from the GRIB message in the Grib Header
- *
- * INTERFACE:
- * int make_grib_log (input_fn, lookup_fn, msg_length, offset,
- * pds, gds, bds, bms, grib_data, errmsg)
- *
- * ARGUMENTS (I=input, O=output, I&O=input and output):
- * (I) char *input_fn; name of input GRIB file
- * (I) char *lookup_fn; name of Lookup file, nil if not used
- * (I) unsigned long msg_length; total length of GRIB message
- * (I) long offset; starting location of GRIB message in bytes
- * (I) PDS_INPUT pds; product definition section structure
- * (I) grid_desc_sec gds; grid description section structure
- * (I) BDS_HEAD_INPUT bds; binary data section header structure
- * (I) BMS_INPUT bms; bit map definition section structure
- * (I) float *grib_data; array of decoded data
- *
- * ACCESSES GLOBAL VARS:
- * int UseTables;
- * set to one if lookup table used
- * CTRS_DEFN db_ctr_tbl[NCTRS];
- * predefined array holding Originating Center info
- * PARM_DEFN db_parm_tbl [MAX_PARM_TBLS * NPARM];
- * predefined arr of Parameter info
- * LVL_DEFN db_lvl_tbl [NLVL];
- * predefined arr of Level info struct
- * MODEL_DEFN db_mdl_tbl [NMODEL];
- * predefined arr of Model info struct
- * GEOM_DEFN db_geom_tbl [NGEOM];
- * predefined arr of Geometry info struct
- *
- * RETURN CODE:
- * 0> no errors; file GRIB.log has been created;
- * 1> error, errmsg filled;
- **********************************************************************
- */
- extern int UseTables; /* set means use lookup tbl defns */
- extern PARM_DEFN db_parm_tbl[]; /* parameter conversion info */
- extern MODEL_DEFN db_mdl_tbl[]; /* model conversion info */
- extern LVL_DEFN db_lvl_tbl[]; /* level conversion info */
- extern GEOM_DEFN db_geom_tbl[]; /* Geom conversion info */
- extern CTR_DEFN db_ctr_tbl[]; /* Ctr conversion info */
- #if PROTOTYPE_NEEDED
- int make_grib_log ( char *input_fn,
- char *lookup_fn,
- unsigned long msg_length,
- long offset,
- PDS_INPUT pds,
- grid_desc_sec gds,
- BDS_HEAD_INPUT bds,
- BMS_INPUT bms,
- float *grib_data,
- char *errmsg)
- #else
- int make_grib_log (input_fn, lookup_fn, msg_length, offset,
- pds, gds, bds, bms, grib_data,errmsg)
- char *input_fn;
- char *lookup_fn;
- unsigned long msg_length;
- long offset;
- PDS_INPUT pds;
- grid_desc_sec gds;
- BDS_HEAD_INPUT bds;
- BMS_INPUT bms;
- float *grib_data;
- char *errmsg;
- #endif
- {
- char *func="make_grib_log";
- int i, indx, k, fd, numpts=100;
- float dsf, res, min, max;
- FILE *fp;
- /*
- *
- * A.0 DEBUG printing
- */
- DPRINT1 ("Entering %s\n", func);
- /*
- *
- * A.1 OPEN file "GRIB.log" in APPEND mode
- */
- fp=fopen ("GRIB.log", "a+");
- if (!fp) {
- DPRINT1("%s: failed to open 'GRIB.log' for appending, skip logfile\n",
- func);
- sprintf (errmsg, "%s: failed to open 'GRIB.log'\n", func);
- return (1);
- }
- /*
- *
- * A.2 WRITE Indicator Section information to file
- * !message length
- * !GRIB Edition number
- */
- fseek(fp, 0L, 2);
- if (ftell(fp) == 0L)
- fprintf (fp, "%s: InFile= %s\n%s: Lookup=%d, fn='%s'\n\n" ,
- func, input_fn, func,UseTables, lookup_fn);
- fprintf (fp, "**** VALID MESSAGE FOUND AT %ld BYTES ****\n" , offset);
- fprintf(fp, "\n********* SECTION 0 IDS *********\n" );
- fprintf(fp, "Total Message length = %ld\nEdition Number = %d\n",
- msg_length, pds.usEd_num);
- /*
- *
- * A.3 WRITE Product Definition Section information to file
- * !Section length
- * !Parameter Table version
- * !Parameter Sub-Table version if defined and flagged by Extension flag
- * !Tracking id if defined and flagged by Extension flag
- */
- fprintf(fp, "\n********* SECTION 1 PDS *********\n" \
- "Section length = %d\nTable version = %d\n",
- pds.uslength, pds.usParm_tbl);
- if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG )
- {
- if (pds.usSub_tbl != 0)
- fprintf(fp,"Local Table version = %d\n",pds.usSub_tbl);
- if(pds.usTrack_num != 0)
- fprintf(fp,"Tracking ID = %d\n",pds.usTrack_num);
- }
-
- /*
- * !Originating Center id
- * !IF (using tables) Name of Originating Center
- */
- fprintf(fp,"Originating Center id = %d\n",pds.usCenter_id);
- if (UseTables)
- if ( db_ctr_tbl[pds.usCenter_id].ctr_dsc[0] )
- fprintf(fp,"Originating Center = %s\n",
- db_ctr_tbl[pds.usCenter_id].ctr_dsc);
- else
- fprintf(fp,"Originating Center ID %d not defined in current table.\n",
- pds.usCenter_id);
- /*
- * !Sub-Table Entry for Originating Center if non-zero and if
- * !extension flag is set
- */
- if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG &&
- pds.usCenter_sub != 0)
- fprintf(fp,"Sub-Table Entry Originating Center = %d\n",pds.usCenter_sub);
- /*
- * !Extension flag
- */
- fprintf(fp,"Extension flag = %d (extensions %s)\n", pds.usExt_flag,
- (pds.usExt_flag == (unsigned short)EXTENSION_FLAG ? "used" : "not used"));
- /*
- * !Model Identification
- * !IF (using tables) Model Description
- */
- fprintf(fp,"Model id = %d\n",pds.usProc_id);
- if (UseTables)
- if ( db_mdl_tbl[pds.usProc_id].grib_dsc[0] )
- fprintf(fp,"Model Description = %s\n",
- db_mdl_tbl[pds.usProc_id].grib_dsc);
- else
- fprintf(fp,"Model ID %d not defined in current table.\n",
- pds.usProc_id);
- /*
- * !Grid Identification
- * !IF (using tables) Grid Description
- */
- fprintf(fp,"Grid id = %d\n",pds.usGrid_id);
- if (UseTables)
- if ( db_geom_tbl[pds.usGrid_id].grib_dsc[0] )
- fprintf(fp,"Grid Description = %s\n",
- db_geom_tbl[pds.usGrid_id].grib_dsc);
- else
- fprintf(fp,"Grid ID %d not defined in current table.\n",
- pds.usGrid_id);
- /*
- * !Parameter Identification
- */
- fprintf(fp,"Parameter id = %d\n",pds.usParm_id);
- /*
- * !IF (usExt_flag is set AND
- * ! (Parm id between 250 and 254) AND (Sub Parm ID defined)))
- * ! PRINT Parm_sub
- * !ENDIF
- */
- if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG &&
- pds.usParm_id>=250 && pds.usParm_id<=254 && pds.usParm_sub!=0)
- fprintf(fp,"Parameter sub-id = %d\n",pds.usParm_sub);
- /*
- * !IF (using lookup table) THEN
- * ! CALCULATE index in Parm Conversion Array to use
- * ! let index= (usParm_Id - 249)*256 + usParm_sub;
- * !
- * ! IF this index in Parm Conversion Array is defined THEN
- * ! PRINT its grib_dsc and grib_unit_dsc
- * ! ELSE
- * ! PRINT it's not defined mesage
- * ! ENDIF
- * !ENDIF
- */
- if(UseTables)
- {
- indx = PARMTBL_INDX (pds.usParm_id, pds.usParm_sub);
- if ( db_parm_tbl[indx].grib_dsc[0] ) {
- fprintf(fp,"Parameter name = %s\n",db_parm_tbl[indx].grib_dsc);
- fprintf(fp,"Parameter units = %s\n",db_parm_tbl[indx].grib_unit_dsc);
- }
- else fprintf(fp,"Parameter ID %d not defined in current table.\n",
- pds.usParm_id);
- }
- /*
- * !Level Id
- * !IF (using tables)
- * ! Level description
- * ! SWITCH (number of octets to store Height1)
- * ! 2: Level = Height1
- * ! 1: Bottom of Layer = Height1
- * ! Top of Layer = Height2
- * ! 0: (no Height value required)
- * ! default: (corrupt table entry or message)
- * ! ENDSWITCH
- * !ELSE (not using tables)
- * ! Level = Height1 (Level assumed)
- * !ENDIF
- */
- fprintf(fp,"Level_type = %d\n",pds.usLevel_id);
- if(UseTables) {
- if ( db_lvl_tbl[pds.usLevel_id].grib_dsc[0] ) {
- fprintf(fp,"Level description = %s\n",
- db_lvl_tbl[pds.usLevel_id].grib_dsc);
- switch(db_lvl_tbl[pds.usLevel_id].num_octets){
- case 2:
- fprintf(fp,"%s = %u\n",
- db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1);
- break;
- case 1:
- fprintf(fp,"%s = %u\n%s = %u\n",
- db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1,
- db_lvl_tbl[pds.usLevel_id].lvl_name_2, pds.usHeight2);
- break;
- case 0:
- break;
- default:
- fprintf(fp,"***Number of octets for table 3 undefined - possibly "
- "corrupt dataset.***\n");
- }
- }else
- fprintf(fp,"Level ID %d not defined in current table.\n",
- pds.usLevel_id);
- } /* end UseTables 'if' statement */
- else fprintf(fp,"Level = %u\n",pds.usHeight1);
- /*
- * !Reference Date/Time:
- * ! Century
- * ! Year
- * ! Month
- * ! Day
- * ! Hour
- * ! Minute
- * ! Second if defined
- */
- fprintf(fp,
- "Reference Date/Time of Data Set:\n" \
- " Century = %d\n Year = %d\n Month = %d\n Day = %d\n"\
- " Hour = %d\n Minute = %d\n",
- pds.usCentury,pds.usYear,pds.usMonth,pds.usDay,pds.usHour,pds.usMinute);
- if(pds.usExt_flag == (unsigned short)EXTENSION_FLAG)
- fprintf(fp," Second = %d\n",pds.usSecond);
- /*
- * !Forecast Time Unit
- * ! Forecast Period 1
- * ! Forecast Period 2
- */
- switch(pds.usFcst_unit_id){
- case 0: fprintf(fp,"Forecast Time Unit = Minute\n"); break;
- case 1: fprintf(fp,"Forecast Time Unit = Hour\n"); break;
- case 2: fprintf(fp,"Forecast Time Unit = Day\n"); break;
- case 3: fprintf(fp,"Forecast Time Unit = Month\n"); break;
- case 4: fprintf(fp,"Forecast Time Unit = Year\n"); break;
- case 5: fprintf(fp,"Forecast Time Unit = Decade (10 years)\n"); break;
- case 6: fprintf(fp,"Forecast Time Unit = Normal (30 years)\n"); break;
- case 7: fprintf(fp,"Forecast Time Unit = Century (100 years)\n"); break;
- case 254: fprintf(fp,"Forecast Time Unit = Second\n"); break;
- default: fprintf(fp,"Forecast Time Unit = UNDEFINED!!\n");
- }
- fprintf(fp," Forecast Period 1 = %d\n",pds.usP1);
- fprintf(fp," Forecast Period 2 = %d\n",pds.usP2);
- /*
- * !Time Range Indicator
- * !Number in Average
- * !Number Missing
- */
- fprintf(fp,"Time Range = %d\n",pds.usTime_range);
- fprintf(fp,"Number in Average = %d\n",pds.usTime_range_avg);
- fprintf(fp,"Number Missing = %d\n",pds.usTime_range_mis);
- /*
- * !Decimal Scale Factor
- */
- fprintf(fp,"Decimal Scale Factor = %d\n",pds.sDec_sc_fctr);
- /*
- *
- * A.4 IF (GDS included) THEN
- * A.4.1 WRITE Grid Definition Section information to file
- * !Section length
- * !Parm_nv
- * !Parm_pv_pl
- * !Data type
- */
- if(pds.usGds_bms_id >> 7 & 1) {
- fprintf(fp,"\n********* SECTION 2 GDS *********\n");
- fprintf(fp,"Section length = %d\n",gds.head.uslength);
- fprintf(fp,"Parm_nv = %d\n",gds.head.usNum_v);
- fprintf(fp,"Parm_pv_pl = %d\n",gds.head.usPl_Pv);
- fprintf(fp,"Data_type = %d\n",gds.head.usData_type);
-
- /*
- * A.4.2 SWITCH (Data Type, Table 6)
- * ! For each Data Type, write the following to file:
- * ! Number of points along rows/columns of grid
- * ! Reference Lat/Lon information
- * ! Resolution and Component Flags (Table 7)
- * ! Direction increments if given
- * ! Assumption of Earth shape
- * ! U&V component orientation
- * ! Scanning mode flags (Table 8)
- * Default: Projection not supported, exit;
- */
- fprintf(fp,"Projection = %s\n", prjn_name[gds.head.usData_type]);
- switch(gds.head.usData_type)
- {
- /*
- * Case 0: Lat/Lon projection
- * Case 10: Rotated Lat/Lon projection
- * Case 20: Stretched Lat/Lon projection
- * Case 30: Stretched Rotated Lat/Lon projection
- */
- case LATLON_PRJ: /* Lat/Lon Grid */
- case ROT_LATLON_PRJ: /* Rotated Lat/Lon */
- case STR_LATLON_PRJ: /* Stretched Lat/Lon */
- case STR_ROT_LATLON_PRJ : /* Stretched and Rotated Lat/Lon */
- fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi);
- fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj);
- fprintf(fp,"Latitude of first grid point = %.3f deg\n",
- ((float)gds.llg.lLat1)/1000.);
- fprintf(fp,"Longitude of first grid point = %.3f deg\n",
- ((float)gds.llg.lLon1)/1000.);
- fprintf(fp,"Latitude of last grid point = %.3f deg\n",
- ((float)gds.llg.lLat2)/1000.);
- fprintf(fp,"Longitude of last grid point = %.3f deg\n",
- ((float)gds.llg.lLon2)/1000.);
- fprintf(fp,"Resolution and Component Flags: \n");
- if ((gds.llg.usRes_flag >> 7) & 1) {
- fprintf(fp," Longitudinal increment = %f deg\n",
- ((float)gds.llg.iDi)/1000.);
- fprintf(fp," Latitudinal increment = %f deg\n",
- ((float)gds.llg.iDj)/1000.);
- }else fprintf(fp," Direction increments not given.\n");
- if ((gds.llg.usRes_flag >> 6) & 1)
- fprintf(fp," Earth assumed oblate spherical.\n");
- else fprintf(fp," Earth assumed spherical.\n");
- if ((gds.llg.usRes_flag >> 3) & 1)
- fprintf(fp," U&V components resolved relative to +I and "
- "+J\n");
- else fprintf(fp," U&V components resolved relative to east "
- "and north.\n");
- fprintf(fp,"Scanning Mode Flags: \n");
- if ((gds.llg.usScan_mode >> 7) & 1)
- fprintf(fp," Points scan in -I direction.\n");
- else fprintf(fp," Points scan in +I direction.\n");
- if ((gds.llg.usScan_mode >> 6) & 1)
- fprintf(fp," Points scan in +J direction.\n");
- else fprintf(fp," Points scan in -J direction.\n");
- if ((gds.llg.usScan_mode >> 5) & 1)
- fprintf(fp," Adjacent points in J direction are "
- "consecutive.\n");
- else fprintf(fp," Adjacent points in I direction are "
- "consecutive.\n");
- /* added 02/98
- This code pertains only to the Stretch/Rotate grids,
- so skip over if it's a LATLON_PRJN type;
- */
- if (gds.head.usData_type != LATLON_PRJ)
- {
- fprintf(fp," Latitude of southern pole = %.3f deg\n",
- ((float)gds.llg.lLat_southpole)/1000.);
- fprintf(fp," Longitude of southern pole = %.3f deg\n",
- ((float)gds.llg.lLon_southpole)/1000.);
- /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
- a single precision floating point value */
- fprintf(fp," Angle of rotation = %.3f\n",
- grib_ibm_local ((unsigned long) gds.llg.lRotate));
- fprintf(fp," Latitude of pole of stretching = %.3f deg\n",
- ((float)gds.llg.lPole_lat)/1000.);
- fprintf(fp," Longitude of pole of stretching = %.3f deg\n",
- ((float)gds.llg.lPole_lon)/1000.);
- /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
- a single precision floating point value */
- fprintf(fp," Stretching factor = %.3f\n",
- grib_ibm_local ((unsigned long)gds.llg.lStretch));
- }
- break;
- /*
- * Case 1: Mercator Projection
- */
- case MERC_PRJ: /* Mercator Projection Grid */
- fprintf(fp,"Number of points along a parallel = %d\n",gds.merc.cols);
- fprintf(fp,"Number of points along a meridian = %d\n",gds.merc.rows);
- fprintf(fp,"Latitude of first grid point = %.3f deg\n",
- ((float)gds.merc.first_lat)/1000.);
- fprintf(fp,"Longitude of first grid point = %.3f deg\n",
- ((float)gds.merc.first_lon)/1000.);
- fprintf(fp,"Latitude of last grid point = %.3f deg\n",
- ((float)gds.merc.La2)/1000.);
- fprintf(fp,"Longitude of last grid point = %.3f deg\n",
- ((float)gds.merc.Lo2)/1000.);
- fprintf(fp,"Latitude of intersection with Earth = %.3f deg\n",
- ((float)gds.merc.latin)/1000.);
- fprintf(fp,"Resolution and Component Flags: \n");
- if ((gds.merc.usRes_flag >> 7) & 1) {
- fprintf(fp," Longitudinal increment = %f deg\n",
- ((float)gds.merc.lon_inc)/1000.);
- fprintf(fp," Latitudinal increment = %f deg\n",
- ((float)gds.merc.lat_inc)/1000.);
- }else fprintf(fp," Direction increments not given.\n");
- if ((gds.merc.usRes_flag >> 6) & 1)
- fprintf(fp," Earth assumed oblate spherical.\n");
- else fprintf(fp," Earth assumed spherical.\n");
- if ((gds.merc.usRes_flag >> 3) & 1)
- fprintf(fp," U&V components resolved relative to +I and "
- "+J\n");
- else fprintf(fp," U&V components resolved relative to east "
- "and north.\n");
- fprintf(fp,"Scanning Mode Flags: \n");
- if ((gds.merc.usScan_mode >> 7) & 1)
- fprintf(fp," Points scan in -I direction.\n");
- else fprintf(fp," Points scan in +I direction.\n");
- if ((gds.merc.usScan_mode >> 6) & 1)
- fprintf(fp," Points scan in +J direction.\n");
- else fprintf(fp," Points scan in -J direction.\n");
- if ((gds.merc.usScan_mode >> 5) & 1)
- fprintf(fp," Adjacent points in J direction are "
- "consecutive.\n");
- else fprintf(fp," Adjacent points in I direction are "
- "consecutive.\n");
- break;
- /*
- * Case 3: Lambert Conformal Projection
- * Case 13: Oblique Lambert Conformal Projection
- * Case 8: Alberts equal-area secant/tangent conic/bipolar Prj
- */
- case LAMB_PRJ: /* Lambert Conformal */
- case OBLIQ_LAMB_PRJ: /* Oblique Lambert Conformal */
- case ALBERS_PRJ: /* Albers equal-area */
-
- fprintf(fp,"Number of points along X-axis = %d\n",gds.lam.iNx);
- fprintf(fp,"Number of points along Y-axis = %d\n",gds.lam.iNy);
- fprintf(fp,"Latitude of first grid point = %.3f deg\n",
- ((float)gds.lam.lLat1)/1000.);
- fprintf(fp,"Longitude of first grid point = %.3f deg\n",
- ((float)gds.lam.lLon1)/1000.);
- fprintf(fp,"Orientation of grid = %.3f deg\n",
- ((float)gds.lam.lLon_orient)/1000.);
- fprintf(fp,"First Latitude Cut = %.3f deg\n",
- ((float)gds.lam.lLat_cut1)/1000.);
- fprintf(fp,"Second Latitude Cut = %.3f deg\n",
- ((float)gds.lam.lLat_cut2)/1000.);
- fprintf(fp,"Resolution and Component Flags: \n");
- if ((gds.lam.usRes_flag >> 7) & 1) {
- fprintf(fp," X-direction increment = %d meters\n",
- gds.lam.ulDx);
- fprintf(fp," Y-direction increment = %d meters\n",
- gds.lam.ulDy);
- }else fprintf(fp," Direction increments not given.\n");
- if ((gds.lam.usRes_flag >> 6) & 1)
- fprintf(fp," Earth assumed oblate spherical.\n");
- else fprintf(fp," Earth assumed spherical.\n");
- if ((gds.lam.usRes_flag >> 3) & 1)
- fprintf(fp," U&V components resolved relative to +I and "
- "+J\n");
- else fprintf(fp," U&V components resolved relative to east "
- "and north.\n");
- fprintf(fp,"Scanning Mode Flags: \n");
- if ((gds.lam.usScan_mode >> 7) & 1)
- fprintf(fp," Points scan in -I direction.\n");
- else fprintf(fp," Points scan in +I direction.\n");
- if ((gds.lam.usScan_mode >> 6) & 1)
- fprintf(fp," Points scan in +J direction.\n");
- else fprintf(fp," Points scan in -J direction.\n");
- if ((gds.lam.usScan_mode >> 5) & 1)
- fprintf(fp," Adjacent points in J direction are "
- "consecutive.\n");
- else fprintf(fp," Adjacent points in I direction are "
- "consecutive.\n");
- /* 02/98 This code pertains only to the Albers projection */
- if (gds.head.usData_type == ALBERS_PRJ)
- {
- fprintf(fp," Latitude of the southern pole = %.3f\n",
- ((float)gds.lam.lLat_southpole)/1000.);
- fprintf(fp," Longitude of the southern pole = %.3f\n",
- ((float)gds.lam.lLon_southpole)/1000.);
- }
- break;
- /*
- * Case 4: Gaussian Lat/Lon Projection
- * Case 14: Rotated Gaussian Lat/Lon Projection
- * Case 24: Stretched Gaussian Lat/Lon Projection
- * Case 34: Stretched Rotated Gaussian Lat/Lon Projection
- */
- case GAUSS_PRJ: /* Gaussian Latitude/Longitude Grid */
- case ROT_GAUSS_PRJ: /* Rotated Gaussian */
- case STR_GAUSS_PRJ : /* Stretched Gaussian */
- case STR_ROT_GAUSS_PRJ : /* Stretched and Rotated Gaussian */
- fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi);
- fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj);
- fprintf(fp,"Latitude of first grid point = %.3f deg\n",
- ((float)gds.llg.lLat1)/1000.);
- fprintf(fp,"Longitude of first grid point = %.3f deg\n",
- ((float)gds.llg.lLon1)/1000.);
- fprintf(fp,"Latitude of last grid point = %.3f deg\n",
- ((float)gds.llg.lLat2)/1000.);
- fprintf(fp,"Longitude of last grid point = %.3f deg\n",
- ((float)gds.llg.lLon2)/1000.);
- fprintf(fp,"Resolution and Component Flags: \n");
- if ((gds.llg.usRes_flag >> 7) & 1) {
- fprintf(fp," i direction increment = %f deg\n",
- ((float)gds.llg.iDi)/1000.);
- fprintf(fp,
- " Number of parallels between pole and equator = %d\n",
- gds.llg.iDj);
- }else fprintf(fp," Direction increments not given.\n");
- if ((gds.llg.usRes_flag >> 6) & 1)
- fprintf(fp," Earth assumed oblate spherical.\n");
- else fprintf(fp," Earth assumed spherical.\n");
- if ((gds.llg.usRes_flag >> 3) & 1)
- fprintf(fp," U&V components resolved relative to +I and "
- "+J\n");
- else fprintf(fp," U&V components resolved relative to east "
- "and north.\n");
- fprintf(fp,"Scanning Mode Flags: \n");
- if ((gds.llg.usScan_mode >> 7) & 1)
- fprintf(fp," Points scan in -I direction.\n");
- else fprintf(fp," Points scan in +I direction.\n");
- if ((gds.llg.usScan_mode >> 6) & 1)
- fprintf(fp," Points scan in +J direction.\n");
- else fprintf(fp," Points scan in -J direction.\n");
- if ((gds.llg.usScan_mode >> 5) & 1)
- fprintf(fp," Adjacent points in J direction are "
- "consecutive.\n");
- else fprintf(fp," Adjacent points in I direction are "
- "consecutive.\n");
- /* added 02/98
- This code pertains only to the Stretch/Rotate grids
- */
- if (gds.head.usData_type != GAUSS_PRJ)
- {
- fprintf(fp," Latitude of southern pole = %.3f deg\n",
- ((float)gds.llg.lLat_southpole)/1000.);
- fprintf(fp," Longitude of southern pole = %.3f deg\n",
- ((float)gds.llg.lLon_southpole)/1000.);
-
- /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
- a single precision floating point value */
- fprintf(fp," Angle of rotation = %.3f\n",
- grib_ibm_local ((unsigned long) gds.llg.lRotate));
- fprintf(fp," Latitude of pole of stretching = %.3f deg\n",
- (float)gds.llg.lPole_lat);
- fprintf(fp," Longitude of pole of stretching = %.3f deg\n",
- (float)gds.llg.lPole_lon);
-
- /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation
- a single precision floating point value */
- fprintf(fp," Stretching factor = %.3f\n",
- grib_ibm_local ((unsigned long)gds.llg.lStretch));
- }
- break;
- /*
- * Case 5: Polar Sterographic Projection
- */
- case POLAR_PRJ: /* Polar Stereographic Projection Grid */
- fprintf(fp,"Number of points along X-axis = %d\n",gds.pol.usNx);
- fprintf(fp,"Number of points along Y-axis = %d\n",gds.pol.usNy);
- fprintf(fp,"Latitude of first grid point = %.3f deg\n",
- ((float)gds.pol.lLat1)/1000.);
- fprintf(fp,"Longitude of first grid point = %.3f deg\n",
- ((float)gds.pol.lLon1)/1000.);
- fprintf(fp,"Orientation of grid = %.3f deg\n",
- ((float)gds.pol.lLon_orient)/1000.);
- fprintf(fp,"Projection Center: ");
- if ((gds.pol.usProj_flag >> 7) & 1)
- fprintf(fp,"South Pole\n");
- else fprintf(fp,"North Pole\n");
- fprintf(fp,"Resolution and Component Flags: \n");
- if ((gds.pol.usRes_flag >> 7) & 1) {
- fprintf(fp," X-direction grid length = %d meters\n",gds.pol.ulDx);
- fprintf(fp," Y-direction grid length = %d meters\n",gds.pol.ulDy);
- }else fprintf(fp," Direction increments not given.\n");
- if ((gds.pol.usRes_flag >> 6) & 1)
- fprintf(fp," Earth assumed oblate spherical.\n");
- else fprintf(fp," Earth assumed spherical.\n");
- if ((gds.pol.usRes_flag >> 3) & 1)
- fprintf(fp," U&V components resolved relative to +I and "
- "+J\n");
- else fprintf(fp," U&V components resolved relative to east "
- "and north.\n");
- fprintf(fp,"Scanning Mode Flags: \n");
- if ((gds.pol.usScan_mode >> 7) & 1)
- fprintf(fp," Points scan in -I direction.\n");
- else fprintf(fp," Points scan in +I direction.\n");
- if ((gds.pol.usScan_mode >> 6) & 1)
- fprintf(fp," Points scan in +J direction.\n");
- else fprintf(fp," Points scan in -J direction.\n");
- if ((gds.pol.usScan_mode >> 5) & 1)
- fprintf(fp," Adjacent points in J direction are "
- "consecutive.\n");
- else fprintf(fp," Adjacent points in I direction are "
- "consecutive.\n");
- break;
- default: /* Bad projection: ignore & continue */
- fprintf(stdout, "\n\n***%s WARNING ***:\nProjection %d is INVALID;\n\n",
- func, gds.head.usData_type);
- fprintf(fp,"================================================\n"\
- "%s: projection %d is not currently implemented\n"\
- "================================================\n",
- func, gds.head.usData_type);
- break;
- /*
- * A.4.2 ENDSWITCH (Data Type)
- */
- } /* Switch */
- } /* gds included */
- /*
- *
- * A.4 ELSE
- * PRINT no Gds message
- * A.4 ENDIF
- */
- else fprintf(fp,"\n******* NO SECTION 2 GDS *********\n" );
- /*
- *
- * A.5 IF (Bitmap Section is present)
- * THEN
- * WRITE Bitmap Section information to file
- * ELSE
- * PRINT no bms mesg
- * ENDIF
- */
- if(pds.usGds_bms_id >> 6 & 1) {
- fprintf(fp,"\n********* SECTION 3 BMS **********\n" );
- fprintf(fp,"Section length = %ld\n", bms.uslength);
- if (bms.uslength <= 6)
- fprintf(fp,"Bitmap is predefined (Not in message).\n");
- else fprintf(fp,"Bitmap is included with message.\n");
- fprintf(fp,"Bitmap ID = %d \n", bms.usBMS_id);
- fprintf(fp,"Number of unused bits = %d\n", bms.usUnused_bits);
- fprintf(fp,"Number of datapoints set = %ld\n", bms.ulbits_set);
- }else{
- fprintf(fp,"\n******* NO SECTION 3 BMS *********\n" );
- }
- /*
- *
- * A.6 WRITE out Binary Data Section Information to file
- * !Section Length
- */
- fprintf(fp,"\n********* SECTION 4 BDS *********\n" );
- fprintf(fp,"Section length = %ld\n",bds.length);
- /*
- * !Table 11 Flags
- */
- fprintf(fp,"Table 11 Flags:\n");
- if ((bds.usBDS_flag >> 7) & 1)
- fprintf(fp," Spherical harmonic coefficients.\n");
- else fprintf(fp," Grid-point data.\n");
- if ((bds.usBDS_flag >> 6) & 1)
- fprintf(fp," Second-order packing.\n");
- else fprintf(fp," Simple Packing.\n");
- if ((bds.usBDS_flag >> 5) & 1)
- fprintf(fp," Integer values.\n");
- else fprintf(fp," Floating point values.\n");
- if ((bds.usBDS_flag >> 4) & 1)
- fprintf(fp," Octet 14 contains additional flag bits.\n");
- else fprintf(fp," No additional flags at octet 14.\n");
- /*
- * !Decimal Scale Factor (Repeated from PDS)
- */
- fprintf(fp,"\nDecimal Scale Factor = %d\n",pds.sDec_sc_fctr);
- /*
- * !Binary Scale Factor
- * !Bit Width
- * !Number of Data Points
- */
- fprintf(fp,"Binary scale factor = %d\n", bds.Bin_sc_fctr);
- fprintf(fp,"Bit width = %d\n", bds.usBit_pack_num);
- fprintf(fp,"Number of data points = %ld\n",bds.ulGrid_size);
- /*
- * A.6.1 WRITE Data Summary to file
- * !Compute Data Min/Max and Resolution
- */
- dsf = (float) pow( (double) 10, (double) pds.sDec_sc_fctr);
- res = (float) pow((double)2,(double)bds.Bin_sc_fctr) / dsf;
- min = bds.fReference / dsf;
- max = (float) (pow((double)2, (double)bds.usBit_pack_num) - 1);
- max = min + max * res;
- fprintf(fp,"Data Minimum = %f\n", min );
- fprintf(fp,"Data Maximum = %f\n", max );
- fprintf(fp,"Resolution = %f\n",res );
- /*
- * !Compute Format Specifier for printing Data
- */
- fd = (int) -1 * (float) log10((double) res) + .5;
- if (fd <= 0)
- {
- fd = 0;
- fprintf(fp,"DATA will be displayed as integers (res > 0.1).\n");
- }
- /*
- * !WRITE First 100 Data Points to file up to Encoded Precision
- */
- if (bds.ulGrid_size > 1) {
- if (bds.ulGrid_size < 100) numpts = bds.ulGrid_size;
- fprintf(fp,"\nDATA ARRAY: (first %d)\n",numpts);
- if (fd > 0) {
- for (i=0; i<numpts; i=i+5)
- if (i + 5 < numpts)
- fprintf(fp, "%03d- %.*f %.*f %.*f %.*f %.*f\n",
- i,fd,grib_data[i],fd,grib_data[i+1],fd,
- grib_data[i+2],fd,grib_data[i+3],fd,grib_data[i+4] );
- else {
- fprintf(fp,"%03d-", i);
- while (i < numpts) fprintf(fp," %.*f", fd, grib_data[i++]);
- fprintf(fp,"\n");
- }
- }
- else {
- for (i=0; i<numpts; i=i+5)
- if (i + 5 < numpts)
- fprintf(fp, "%03d- %.0f %.0f %.0f %.0f %.0f\n",
- i,grib_data[i],grib_data[i+1],
- grib_data[i+2],grib_data[i+3],grib_data[i+4] );
- else {
- fprintf(fp,"%03d-", i);
- while (i < numpts) fprintf(fp," %.0f", grib_data[i++]);
- fprintf(fp,"\n");
- }
- }
-
- }
- fprintf (fp,"\n******** END OF MESSAGE *******\n\n");
- /*
- *
- * A.7 CLOSE file
- */
- fclose (fp);
- /*
- *
- * A.8 DEBUG printing
- */
- DPRINT1 ("Leaving %s, no errors\n", func);
- return (0);
- /*
- * END OF FUNCTION
- *
- *
- */
- }