PageRenderTime 46ms CodeModel.GetById 14ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/external/io_grib1/MEL_grib1/make_grib_log.c

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