Context Navigation

make_grib_log.c

Last change on this file was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 30.1 KB

Line
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	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/WRF.COMMON/WRFV2/external/io_grib1/MEL_grib1/make_grib_log.c

Download in other formats: