?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
updated top-level README and version_decl for V4.4.2 (#1795)
[WRF.git] / external / io_grib1 / WGRIB / wgrib_main.c
blob1f9fbe7ffe07186eca53b47f1e69dc5e1a8c3f72
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <stddef.h>
5 #include <math.h>
6 #include <float.h>
7
8 #include "pds4.h"
9 #include "gds.h"
10 #include "bms.h"
11 #include "bds.h"
12 #include "cnames.h"
13 #include "grib.h"
14
15 #define VERSION "v1.8.0.9i (12-03-04) Wesley Ebisuzaki\n\t\tDWD-tables 2,201-203 (8-19-2003) Helmut P. Frank\n\t\tspectral: Luis Kornblueh (MPI)"
16
17 #define CHECK_GRIB
18
19 /*
20 * wgrib.c extract/inventory grib records
21 *
22 * Wesley Ebisuzaki
23 *
24 * 11/94 - v1.0
25 * 11/94 - v1.1: arbitary size grids, -i option
26 * 11/94 - v1.2: bug fixes, ieee option, more info
27 * 1/95 - v1.2.4: fix headers for SUN acc
28 * 2/95 - v1.2.5: add num_ave in -s listing
29 * 2/95 - v1.2.6: change %d to %ld
30 * 2/95 - v1.2.7: more output, added some polar stereographic support
31 * 2/95 - v1.2.8: max min format changed %f to %g, tidying up more info
32 * 3/95 - v1.3.0: fix bug with bitmap, allow numbers > UNDEFINED
33 * 3/95 - v1.3.1: print number of missing points (verbose)
34 * 3/95 - v1.3.2: -append option added
35 * 4/95 - v1.3.2a,b: more output, polar stereo support (-V option)
36 * 4/95 - v1.3.3: added ECMWF parameter table (prelim)
37 * 6/95 - v1.3.4: nxny from BDS rather than gds?
38 * 9/95 - v1.3.4d: speedup in grib write
39 * 11/95 - v1.3.4f: new ECMWF parameter table (from Mike Fiorino), EC logic
40 * 2/96 - v1.3.4g-h: prelim fix for GDS-less grib files
41 * 2/96 - v1.3.4i: faster missing(), -V: "pos n" -> "n" (field 2)
42 * 3/96 - v1.4: fix return code (!inventory), and short records near EOF
43 * 6/96 - v1.4.1a: faster grib->binary decode, updated ncep parameter table, mod. in clim. desc
44 * 7/96 - v1.5.0: parameter-table aware, -v option changed, added "comments"
45 * increased NTRY to 100 in seek_grib
46 * 11/96 - v1.5.0b: added ECMWF parameter table 128
47 * 1/97 - v1.5.0b2: if nxny != nx*ny { nx = nxny; ny = 1 }
48 * 3/97 - v1.5.0b5: added: -PDS -GDS, Lambert Conformal
49 * 3/97 - v1.5.0b6: added: -verf
50 * 4/97 - v1.5.0b7: added -PDS10, -GDS10 and enhanced -PDS -GDS
51 * 4/97 - v1.5.0b8: "bitmap missing x" -> "bitmap: x undef"
52 * 5/97 - v1.5.0b9: thinned grids meta data
53 * 5/97 - v1.5.0b10: changed 0hr fcst to anal for TR=10 and P1=P2=0
54 * 5/97 - v1.5.0b10: added -H option
55 * 6/97 - v1.5.0b12: thinned lat-long grids -V option
56 * 6/97 - v1.5.0b13: -4yr
57 * 6/97 - v1.5.0b14: fix century mark Y=100 not 0
58 * 7/97 - v1.5.0b15: add ncep opn grib table
59 * 12/97 - v1.6.1.a: made ncep_opn the default table
60 * 12/97 - v1.6.1.b: changed 03TOT to O3TOT in operational ncep table
61 * 1/98 - v1.6.2: added Arakawa E grid meta-data
62 * 1/98 - v1.6.2.1: added some mode data, Scan -> scan
63 * 4/98 - v1.6.2.4: reanalysis id code: subcenter==0 && process==180
64 * 5/98 - v1.6.2.5: fix -H code to write all of GDS
65 * 7/98 - v1.7: fix decoding bug for bitmap and no. bits > 24 (theoretical bug)
66 * 7/98 - v1.7.0.b1: add km to Mercator meta-data
67 * 5/99 - v1.7.2: bug with thinned grids & bitmaps (nxny != nx*ny)
68 * 5/99 - v1.7.3: updated NCEP opn grib table
69 * 8/99 - v1.7.3.1: updated level information
70 * 9/00 - v1.7.3.4a: check for missing grib file
71 * 2/01 - v1.7.3.5: handle data with precision greater than 31 bits
72 * 8/01 - vDWD : added DWD GRIB tables 201, 202, 203, Helmut P. Frank
73 * 9/01 - vDWD : added output "Triangular grid", Helmut P. Frank
74 * 9/01 - v1.7.4: merged Hemut P. Frank's changes to current wgrib source code
75 * 3/02 - vMPIfM: added support for spectral data type
76 * 4/02 - v1.8: merge vMPIfM changes, some fixes/generalizations
77 * 10/02 - v1.8.0.1: added cptec table 254
78 * 10/02 - v1.8.0.2: no test of grib test if no gds, level 117 redone
79 * 10/02 - v1.8.0.3: update ncep_opn grib and levels
80 * 11/02 - v1.8.0.3a: updated ncep_opn and ncep table 129
81 * 9/03 - v1.8.0.4: update dwd tables (Helmut P. Frank), -dwdgrib option
82 * 9/03 - v1.8.0.5: fix scan mode and change format
83 * 10/03 - v1.8.0.7: Changes from Norwegian Met. Inst (ec tab #131, ex_ext)
84 * 10/03 - v1.8.0.8: added -ncep_ens option
85 *
86 */
87
88 /*
89 * MSEEK = I/O buffer size for seek_grib
90 */
91
92 #define MSEEK 1024
93 #define BUFF_ALLOC0 40000
94
95
96 #ifndef min
97 #define min(a,b) ((a) < (b) ? (a) : (b))
98 #define max(a,b) ((a) < (b) ? (b) : (a))
99 #endif
100
101 #ifndef DEF_T62_NCEP_TABLE
102 #define DEF_T62_NCEP_TABLE rean
103 #endif
104 enum Def_NCEP_Table def_ncep_table = DEF_T62_NCEP_TABLE;
105 int minute = 0;
106 int ncep_ens = 0;
107
108 int main(int argc, char **argv) {
109
110 unsigned char *buffer;
111 float *array;
112 double temp, rmin, rmax;
113 int i, nx, ny, file_arg;
114 long int len_grib, pos = 0, nxny, buffer_size, n_dump, count = 1;
115 unsigned char *msg, *pds, *gds, *bms, *bds, *pointer;
116 FILE *input, *dump_file = NULL;
117 char line[200];
118 enum {BINARY, TEXT, IEEE, GRIB, NONE} output_type = NONE;
119 enum {DUMP_ALL, DUMP_RECORD, DUMP_POSITION, DUMP_LIST, INVENTORY}
120 mode = INVENTORY;
121 enum {none, dwd, simple} header = simple;
122
123 long int dump = -1;
124 int verbose = 0, append = 0, v_time = 0, year_4 = 0, output_PDS_GDS = 0;
125 int print_GDS = 0, print_GDS10 = 0, print_PDS = 0, print_PDS10 = 0;
126 char *dump_file_name = "dump", open_parm[3];
127 int return_code = 0;
128
129 if (argc == 1) {
130 fprintf(stderr, "\nPortable Grib decoder for %s etc.\n",
131 (def_ncep_table == opn_nowarn || def_ncep_table == opn) ?
132 "NCEP Operations" : "NCEP/NCAR Reanalysis");
133 fprintf(stderr, " it slices, dices %s\n", VERSION);
134 fprintf(stderr, " usage: %s [grib file] [options]\n\n", argv[0]);
135
136 fprintf(stderr, "Inventory/diagnostic-output selections\n");
137 fprintf(stderr, " -s/-v short/verbose inventory\n");
138 fprintf(stderr, " -V diagnostic output (not inventory)\n");
139 fprintf(stderr, " (none) regular inventory\n");
140
141 fprintf(stderr, " Options\n");
142 fprintf(stderr, " -PDS/-PDS10 print PDS in hex/decimal\n");
143 fprintf(stderr, " -GDS/-GDS10 print GDS in hex/decimal\n");
144 fprintf(stderr, " -verf print forecast verification time\n");
145 fprintf(stderr, " -ncep_opn/-ncep_rean default T62 NCEP grib table\n");
146 fprintf(stderr, " -4yr print year using 4 digits\n");
147 fprintf(stderr, " -min print minutes\n");
148 fprintf(stderr, " -ncep_ens ensemble info encoded in ncep format\n");
149
150 fprintf(stderr, "Decoding GRIB selection\n");
151 fprintf(stderr, " -d [record number|all] decode record number\n");
152 fprintf(stderr, " -p [byte position] decode record at byte position\n");
153 fprintf(stderr, " -i decode controlled by stdin (inventory list)\n");
154 fprintf(stderr, " (none) no decoding\n");
155
156 fprintf(stderr, " Options\n");
157 fprintf(stderr, " -text/-ieee/-grib/-bin convert to text/ieee/grib/bin (default)\n");
158 fprintf(stderr, " -nh/-h output will have no headers/headers (default)\n");
159 fprintf(stderr, " -dwdgrib output dwd headers, grib (do not append)\n");
160 fprintf(stderr, " -H output will include PDS and GDS (-bin/-ieee only)\n");
161 fprintf(stderr, " -append append to output file\n");
162 fprintf(stderr, " -o [file] output file name, 'dump' is default\n");
163 exit(8);
164 }
165 file_arg = 0;
166 for (i = 1; i < argc; i++) {
167 if (strcmp(argv[i],"-PDS") == 0) {
168 print_PDS = 1;
169 continue;
170 }
171 if (strcmp(argv[i],"-PDS10") == 0) {
172 print_PDS10 = 1;
173 continue;
174 }
175 if (strcmp(argv[i],"-GDS") == 0) {
176 print_GDS = 1;
177 continue;
178 }
179 if (strcmp(argv[i],"-GDS10") == 0) {
180 print_GDS10 = 1;
181 continue;
182 }
183 if (strcmp(argv[i],"-v") == 0) {
184 verbose = 1;
185 continue;
186 }
187 if (strcmp(argv[i],"-V") == 0) {
188 verbose = 2;
189 continue;
190 }
191 if (strcmp(argv[i],"-s") == 0) {
192 verbose = -1;
193 continue;
194 }
195 if (strcmp(argv[i],"-text") == 0) {
196 output_type = TEXT;
197 continue;
198 }
199 if (strcmp(argv[i],"-bin") == 0) {
200 output_type = BINARY;
201 continue;
202 }
203 if (strcmp(argv[i],"-ieee") == 0) {
204 output_type = IEEE;
205 continue;
206 }
207 if (strcmp(argv[i],"-grib") == 0) {
208 output_type = GRIB;
209 continue;
210 }
211 if (strcmp(argv[i],"-nh") == 0) {
212 header = none;
213 continue;
214 }
215 if (strcmp(argv[i],"-h") == 0) {
216 header = simple;
217 continue;
218 }
219 if (strcmp(argv[i],"-dwdgrib") == 0) {
220 header = dwd;
221 output_type = GRIB;
222 continue;
223 }
224 if (strcmp(argv[i],"-append") == 0) {
225 append = 1;
226 continue;
227 }
228 if (strcmp(argv[i],"-verf") == 0) {
229 v_time = 1;
230 continue;
231 }
232 if (strcmp(argv[i],"-d") == 0) {
233 if (strcmp(argv[i+1],"all") == 0) {
234 mode = DUMP_ALL;
235 }
236 else {
237 dump = atol(argv[i+1]);
238 mode = DUMP_RECORD;
239 }
240 i++;
241 if (output_type == NONE) output_type = BINARY;
242 continue;
243 }
244 if (strcmp(argv[i],"-p") == 0) {
245 pos = atol(argv[i+1]);
246 i++;
247 dump = 1;
248 if (output_type == NONE) output_type = BINARY;
249 mode = DUMP_POSITION;
250 continue;
251 }
252 if (strcmp(argv[i],"-i") == 0) {
253 if (output_type == NONE) output_type = BINARY;
254 mode = DUMP_LIST;
255 continue;
256 }
257 if (strcmp(argv[i],"-H") == 0) {
258 output_PDS_GDS = 1;
259 continue;
260 }
261 if (strcmp(argv[i],"-NH") == 0) {
262 output_PDS_GDS = 0;
263 continue;
264 }
265 if (strcmp(argv[i],"-4yr") == 0) {
266 year_4 = 1;
267 continue;
268 }
269 if (strcmp(argv[i],"-ncep_opn") == 0) {
270 def_ncep_table = opn_nowarn;
271 continue;
272 }
273 if (strcmp(argv[i],"-ncep_rean") == 0) {
274 def_ncep_table = rean_nowarn;
275 continue;
276 }
277 if (strcmp(argv[i],"-o") == 0) {
278 dump_file_name = argv[i+1];
279 i++;
280 continue;
281 }
282 if (strcmp(argv[i],"--v") == 0) {
283 printf("wgrib: %s\n", VERSION);
284 exit(0);
285 }
286 if (strcmp(argv[i],"-min") == 0) {
287 minute = 1;
288 continue;
289 }
290 if (strcmp(argv[i],"-ncep_ens") == 0) {
291 ncep_ens = 1;
292 continue;
293 }
294 if (file_arg == 0) {
295 file_arg = i;
296 }
297 else {
298 fprintf(stderr,"argument: %s ????\n", argv[i]);
299 }
300 }
301 if (file_arg == 0) {
302 fprintf(stderr,"no GRIB file to process\n");
303 exit(8);
304 }
305 if ((input = fopen(argv[file_arg],"rb")) == NULL) {
306 fprintf(stderr,"could not open file: %s\n", argv[file_arg]);
307 exit(7);
308 }
309
310 if ((buffer = (unsigned char *) malloc(BUFF_ALLOC0)) == NULL) {
311 fprintf(stderr,"not enough memory\n");
312 }
313 buffer_size = BUFF_ALLOC0;
314
315 /* open output file */
316 if (mode != INVENTORY) {
317 open_parm[0] = append ? 'a' : 'w'; open_parm[1] = 'b'; open_parm[2] = '\0';
318 if (output_type == TEXT) open_parm[1] = '\0';
319
320 if ((dump_file = fopen(dump_file_name,open_parm)) == NULL) {
321 fprintf(stderr,"could not open dump file\n");
322 exit(8);
323 }
324 if (header == dwd && output_type == GRIB) wrtieee_header(0, dump_file);
325 }
326
327 /* skip dump - 1 records */
328 for (i = 1; i < dump; i++) {
329 msg = seek_grib(input, &pos, &len_grib, buffer, MSEEK);
330 if (msg == NULL) {
331 fprintf(stderr, "ran out of data or bad file\n");
332 exit(8);
333 }
334 pos += len_grib;
335 }
336 if (dump > 0) count += dump - 1;
337 n_dump = 0;
338
339 for (;;) {
340 if (n_dump == 1 && (mode == DUMP_RECORD || mode == DUMP_POSITION)) break;
341 if (mode == DUMP_LIST) {
342 if (fgets(line,sizeof(line), stdin) == NULL) break;
343 line[sizeof(line) - 1] = 0;
344 if (sscanf(line,"%ld:%ld:", &count, &pos) != 2) {
345 fprintf(stderr,"bad input from stdin\n");
346 fprintf(stderr," %s\n", line);
347 exit(8);
348 }
349 }
350
351 msg = seek_grib(input, &pos, &len_grib, buffer, MSEEK);
352 if (msg == NULL) {
353 if (mode == INVENTORY || mode == DUMP_ALL) break;
354 fprintf(stderr,"missing GRIB record(s)\n");
355 exit(8);
356 }
357
358 /* read all whole grib record */
359 if (len_grib + msg - buffer > buffer_size) {
360 buffer_size = len_grib + msg - buffer + 1000;
361 buffer = (unsigned char *) realloc((void *) buffer, buffer_size);
362 if (buffer == NULL) {
363 fprintf(stderr,"ran out of memory\n");
364 exit(8);
365 }
366 }
367 read_grib(input, pos, len_grib, buffer);
368
369 /* parse grib message */
370
371 msg = buffer;
372 pds = (msg + 8);
373 pointer = pds + PDS_LEN(pds);
374 #ifdef DEBUG
375 printf("LEN_GRIB= 0x%x\n", len_grib);
376 printf("PDS_LEN= 0x%x: at 0x%x\n", PDS_LEN(pds),pds-msg);
377 #endif
378 if (PDS_HAS_GDS(pds)) {
379 gds = pointer;
380 pointer += GDS_LEN(gds);
381 #ifdef DEBUG
382 printf("GDS_LEN= 0x%x: at 0x%x\n", GDS_LEN(gds), gds-msg);
383 #endif
384 }
385 else {
386 gds = NULL;
387 }
388
389 if (PDS_HAS_BMS(pds)) {
390 bms = pointer;
391 pointer += BMS_LEN(bms);
392 #ifdef DEBUG
393 printf("BMS_LEN= 0x%x: at 0x%x\n", BMS_LEN(bms),bms-msg);
394 #endif
395 }
396 else {
397 bms = NULL;
398 }
399
400 bds = pointer;
401 pointer += BDS_LEN(bds);
402 #ifdef DEBUG
403 printf("BDS_LEN= 0x%x: at 0x%x\n", BDS_LEN(bds),bds-msg);
404 #endif
405
406 #ifdef DEBUG
407 printf("END_LEN= 0x%x: at 0x%x\n", 4,pointer-msg);
408 if (pointer-msg+4 != len_grib) {
409 fprintf(stderr,"Len of grib message is inconsistent.\n");
410 }
411 #endif
412
413 /* end section - "7777" in ascii */
414 if (pointer[0] != 0x37 || pointer[1] != 0x37 ||
415 pointer[2] != 0x37 || pointer[3] != 0x37) {
416 fprintf(stderr,"\n\n missing end section\n");
417 fprintf(stderr, "%2x %2x %2x %2x\n", pointer[0], pointer[1],
418 pointer[2], pointer[3]);
419 #ifdef DEBUG
420 printf("ignoring missing end section\n");
421 #else
422 exit(8);
423 #endif
424 }
425
426 /* figure out size of array */
427 if (gds != NULL) {
428 GDS_grid(gds, bds, &nx, &ny, &nxny);
429 }
430 else if (bms != NULL) {
431 nxny = nx = BMS_nxny(bms);
432 ny = 1;
433 }
434 else {
435 if (BDS_NumBits(bds) == 0) {
436 nxny = nx = 1;
437 fprintf(stderr,"Missing GDS, constant record .. cannot "
438 "determine number of data points\n");
439 }
440 else {
441 nxny = nx = BDS_NValues(bds);
442 }
443 ny = 1;
444 }
445
446 #ifdef CHECK_GRIB
447 if (gds && ! GDS_Harmonic(gds)) {
448 /* this grib check only works for simple packing */
449 /* turn off if harmonic */
450 if (BDS_NumBits(bds) != 0) {
451 i = BDS_NValues(bds);
452 if (bms != NULL) {
453 i += missing_points(BMS_bitmap(bms),nxny);
454 }
455 if (i != nxny) {
456 fprintf(stderr,"grib header at record %ld: two values of nxny %ld %d\n",
457 count,nxny,i);
458 fprintf(stderr," LEN %d DataStart %d UnusedBits %d #Bits %d nxny %ld\n",
459 BDS_LEN(bds), BDS_DataStart(bds),BDS_UnusedBits(bds),
460 BDS_NumBits(bds), nxny);
461 return_code = 15;
462 nxny = nx = i;
463 ny = 1;
464 }
465 }
466
467 }
468 #endif
469
470 if (verbose <= 0) {
471 printf("%ld:%ld:d=", count, pos);
472 PDS_date(pds,year_4,v_time);
473 printf(":%s:", k5toa(pds));
474
475 if (verbose == 0) printf("kpds5=%d:kpds6=%d:kpds7=%d:TR=%d:P1=%d:P2=%d:TimeU=%d:",
476 PDS_PARAM(pds),PDS_KPDS6(pds),PDS_KPDS7(pds),
477 PDS_TimeRange(pds),PDS_P1(pds),PDS_P2(pds),
478 PDS_ForecastTimeUnit(pds));
479 levels(PDS_KPDS6(pds), PDS_KPDS7(pds),PDS_Center(pds)); printf(":");
480 PDStimes(PDS_TimeRange(pds),PDS_P1(pds),PDS_P2(pds),
481 PDS_ForecastTimeUnit(pds));
482 if (PDS_Center(pds) == ECMWF) EC_ext(pds,"",":");
483 ensemble(pds, verbose);
484 printf("NAve=%d",PDS_NumAve(pds));
485 if (print_PDS || print_PDS10) print_pds(pds, print_PDS, print_PDS10, verbose);
486 if (gds && (print_GDS || print_GDS10)) print_gds(gds, print_GDS, print_GDS10, verbose);
487 printf("\n");
488 }
489 else if (verbose == 1) {
490 printf("%ld:%ld:D=", count, pos);
491 PDS_date(pds, 1, v_time);
492 printf(":%s:", k5toa(pds));
493 levels(PDS_KPDS6(pds), PDS_KPDS7(pds), PDS_Center(pds)); printf(":");
494 printf("kpds=%d,%d,%d:",
495 PDS_PARAM(pds),PDS_KPDS6(pds),PDS_KPDS7(pds));
496 PDStimes(PDS_TimeRange(pds),PDS_P1(pds),PDS_P2(pds),
497 PDS_ForecastTimeUnit(pds));
498 if (PDS_Center(pds) == ECMWF) EC_ext(pds,"",":");
499 ensemble(pds, verbose);
500 GDS_winds(gds, verbose);
501 printf("\"%s", k5_comments(pds));
502 if (print_PDS || print_PDS10) print_pds(pds, print_PDS, print_PDS10, verbose);
503 if (gds && (print_GDS || print_GDS10)) print_gds(gds, print_GDS, print_GDS10, verbose);
504 printf("\n");
505 }
506 else if (verbose == 2) {
507 printf("rec %ld:%ld:date ", count, pos);
508 PDS_date(pds, 1, v_time);
509 printf(" %s kpds5=%d kpds6=%d kpds7=%d levels=(%d,%d) grid=%d ",
510 k5toa(pds), PDS_PARAM(pds), PDS_KPDS6(pds), PDS_KPDS7(pds),
511 PDS_LEVEL1(pds), PDS_LEVEL2(pds), PDS_Grid(pds));
512 levels(PDS_KPDS6(pds),PDS_KPDS7(pds),PDS_Center(pds));
513
514 printf(" ");
515 if (PDS_Center(pds) == ECMWF) EC_ext(pds,""," ");
516 ensemble(pds, verbose);
517 PDStimes(PDS_TimeRange(pds),PDS_P1(pds),PDS_P2(pds),
518 PDS_ForecastTimeUnit(pds));
519 if (bms != NULL)
520 printf(" bitmap: %d undef", missing_points(BMS_bitmap(bms),nxny));
521 printf("\n %s=%s\n", k5toa(pds), k5_comments(pds));
522
523 printf(" timerange %d P1 %d P2 %d TimeU %d nx %d ny %d GDS grid %d "
524 "num_in_ave %d missing %d\n",
525 PDS_TimeRange(pds),PDS_P1(pds),PDS_P2(pds),
526 PDS_ForecastTimeUnit(pds), nx, ny,
527 gds == NULL ? -1 : GDS_DataType(gds),
528 PDS_NumAve(pds), PDS_NumMissing(pds));
529
530 printf(" center %d subcenter %d process %d Table %d",
531 PDS_Center(pds),PDS_Subcenter(pds),PDS_Model(pds),
532 PDS_Vsn(pds));
533 GDS_winds(gds, verbose);
534 printf("\n");
535
536 if (gds && GDS_LatLon(gds) && nx != -1)
537 printf(" latlon: lat %f to %f by %f nxny %ld\n"
538 " long %f to %f by %f, (%d x %d) scan %d "
539 "mode %d bdsgrid %d\n",
540 0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
541 0.001*GDS_LatLon_dy(gds), nxny, 0.001*GDS_LatLon_Lo1(gds),
542 0.001*GDS_LatLon_Lo2(gds), 0.001*GDS_LatLon_dx(gds),
543 nx, ny, GDS_LatLon_scan(gds), GDS_LatLon_mode(gds),
544 BDS_Grid(bds));
545 else if (gds && GDS_LatLon(gds) && nx == -1) {
546 printf(" thinned latlon: lat %f to %f by %f nxny %ld\n"
547 " long %f to %f, %ld grid pts (%d x %d) scan %d"
548 " mode %d bdsgrid %d\n",
549 0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
550 0.001*GDS_LatLon_dy(gds), nxny, 0.001*GDS_LatLon_Lo1(gds),
551 0.001*GDS_LatLon_Lo2(gds),
552 nxny, nx, ny, GDS_LatLon_scan(gds), GDS_LatLon_mode(gds),
553 BDS_Grid(bds));
554 GDS_prt_thin_lon(gds);
555 }
556 else if (gds && GDS_Gaussian(gds) && nx != -1)
557 printf(" gaussian: lat %f to %f\n"
558 " long %f to %f by %f, (%d x %d) scan %d"
559 " mode %d bdsgrid %d\n",
560 0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
561 0.001*GDS_LatLon_Lo1(gds), 0.001*GDS_LatLon_Lo2(gds),
562 0.001*GDS_LatLon_dx(gds),
563 nx, ny, GDS_LatLon_scan(gds), GDS_LatLon_mode(gds),
564 BDS_Grid(bds));
565 else if (gds && GDS_Gaussian(gds) && nx == -1) {
566 printf(" thinned gaussian: lat %f to %f\n"
567 " lon %f %ld grid pts (%d x %d) scan %d"
568 " mode %d bdsgrid %d nlat:\n",
569 0.001*GDS_LatLon_La1(gds), 0.001*GDS_LatLon_La2(gds),
570 0.001*GDS_LatLon_Lo1(gds),
571 nxny, nx, ny, GDS_LatLon_scan(gds), GDS_LatLon_mode(gds),
572 BDS_Grid(bds));
573 GDS_prt_thin_lon(gds);
574 }
575 else if (gds && GDS_Polar(gds))
576 printf(" polar stereo: Lat1 %f Long1 %f Orient %f\n"
577 " %s pole (%d x %d) Dx %d Dy %d scan %d mode %d\n",
578 0.001*GDS_Polar_La1(gds),0.001*GDS_Polar_Lo1(gds),
579 0.001*GDS_Polar_Lov(gds),
580 GDS_Polar_pole(gds) == 0 ? "north" : "south", nx,ny,
581 GDS_Polar_Dx(gds),GDS_Polar_Dy(gds),
582 GDS_Polar_scan(gds), GDS_Polar_mode(gds));
583 else if (gds && GDS_Lambert(gds))
584 printf(" Lambert Conf: Lat1 %f Lon1 %f Lov %f\n"
585 " Latin1 %f Latin2 %f LatSP %f LonSP %f\n"
586 " %s (%d x %d) Dx %f Dy %f scan %d mode %d\n",
587 0.001*GDS_Lambert_La1(gds),0.001*GDS_Lambert_Lo1(gds),
588 0.001*GDS_Lambert_Lov(gds),
589 0.001*GDS_Lambert_Latin1(gds), 0.001*GDS_Lambert_Latin2(gds),
590 0.001*GDS_Lambert_LatSP(gds), 0.001*GDS_Lambert_LonSP(gds),
591 GDS_Lambert_NP(gds) ? "North Pole": "South Pole",
592 GDS_Lambert_nx(gds), GDS_Lambert_ny(gds),
593 0.001*GDS_Lambert_dx(gds), 0.001*GDS_Lambert_dy(gds),
594 GDS_Lambert_scan(gds), GDS_Lambert_mode(gds));
595 else if (gds && GDS_Mercator(gds))
596 printf(" Mercator: lat %f to %f by %f km nxny %ld\n"
597 " long %f to %f by %f km, (%d x %d) scan %d"
598 " mode %d Latin %f bdsgrid %d\n",
599 0.001*GDS_Merc_La1(gds), 0.001*GDS_Merc_La2(gds),
600 0.001*GDS_Merc_dy(gds), nxny, 0.001*GDS_Merc_Lo1(gds),
601 0.001*GDS_Merc_Lo2(gds), 0.001*GDS_Merc_dx(gds),
602 nx, ny, GDS_Merc_scan(gds), GDS_Merc_mode(gds),
603 0.001*GDS_Merc_Latin(gds), BDS_Grid(bds));
604 else if (gds && GDS_ssEgrid(gds))
605 printf(" Semi-staggered Arakawa E-Grid: lat0 %f lon0 %f nxny %d\n"
606 " dLat %f dLon %f (%d x %d) scan %d mode %d\n",
607 0.001*GDS_ssEgrid_La1(gds), 0.001*GDS_ssEgrid_Lo1(gds),
608 GDS_ssEgrid_n(gds)*GDS_ssEgrid_n_dum(gds),
609 0.001*GDS_ssEgrid_dj(gds), 0.001*GDS_ssEgrid_di(gds),
610 GDS_ssEgrid_Lo2(gds), GDS_ssEgrid_La2(gds),
611 GDS_ssEgrid_scan(gds), GDS_ssEgrid_mode(gds));
612 else if (gds && GDS_ss2dEgrid(gds))
613 printf(" Semi-staggered Arakawa E-Grid (2D): lat0 %f lon0 %f nxny %d\n"
614 " dLat %f dLon %f (tlm0d %f tph0d %f) scan %d mode %d\n",
615 0.001*GDS_ss2dEgrid_La1(gds), 0.001*GDS_ss2dEgrid_Lo1(gds),
616 GDS_ss2dEgrid_nx(gds)*GDS_ss2dEgrid_ny(gds),
617 0.001*GDS_ss2dEgrid_dj(gds), 0.001*GDS_ss2dEgrid_di(gds),
618 0.001*GDS_ss2dEgrid_Lo2(gds), 0.001*GDS_ss2dEgrid_La2(gds),
619 GDS_ss2dEgrid_scan(gds), GDS_ss2dEgrid_mode(gds));
620 else if (gds && GDS_fEgrid(gds))
621 printf(" filled Arakawa E-Grid: lat0 %f lon0 %f nxny %d\n"
622 " dLat %f dLon %f (%d x %d) scan %d mode %d\n",
623 0.001*GDS_fEgrid_La1(gds), 0.001*GDS_fEgrid_Lo1(gds),
624 GDS_fEgrid_n(gds)*GDS_fEgrid_n_dum(gds),
625 0.001*GDS_fEgrid_dj(gds), 0.001*GDS_fEgrid_di(gds),
626 GDS_fEgrid_Lo2(gds), GDS_fEgrid_La2(gds),
627 GDS_fEgrid_scan(gds), GDS_fEgrid_mode(gds));
628 else if (gds && GDS_RotLL(gds))
629 printf(" rotated LatLon grid lat %f to %f lon %f to %f\n"
630 " nxny %ld (%d x %d) dx %d dy %d scan %d mode %d\n"
631 " transform: south pole lat %f lon %f rot angle %f\n",
632 0.001*GDS_RotLL_La1(gds), 0.001*GDS_RotLL_La2(gds),
633 0.001*GDS_RotLL_Lo1(gds), 0.001*GDS_RotLL_Lo2(gds),
634 nxny, GDS_RotLL_nx(gds), GDS_RotLL_ny(gds),
635 GDS_RotLL_dx(gds), GDS_RotLL_dy(gds),
636 GDS_RotLL_scan(gds), GDS_RotLL_mode(gds),
637 0.001*GDS_RotLL_LaSP(gds), 0.001*GDS_RotLL_LoSP(gds),
638 GDS_RotLL_RotAng(gds) );
639 else if (gds && GDS_Gnomonic(gds))
640 printf(" Gnomonic grid\n");
641 else if (gds && GDS_Harmonic(gds))
642 printf(" Harmonic (spectral): pentagonal spectral truncation: nj %d nk %d nm %d\n",
643 GDS_Harmonic_nj(gds), GDS_Harmonic_nk(gds),
644 GDS_Harmonic_nm(gds));
645 if (gds && GDS_Harmonic_type(gds) == 1)
646 printf(" Associated Legendre polynomials\n");
647 else if (gds && GDS_Triangular(gds))
648 printf(" Triangular grid: nd %d ni %d (= 2^%d x 3^%d)\n",
649 GDS_Triangular_nd(gds), GDS_Triangular_ni(gds),
650 GDS_Triangular_ni2(gds), GDS_Triangular_ni3(gds) );
651 if (print_PDS || print_PDS10)
652 print_pds(pds, print_PDS, print_PDS10, verbose);
653 if (gds && (print_GDS || print_GDS10))
654 print_gds(gds, print_GDS, print_GDS10, verbose);
655 }
656
657 if (mode != INVENTORY && output_type == GRIB) {
658 if (header == dwd) wrtieee_header((int) len_grib, dump_file);
659 fwrite((void *) msg, sizeof(char), len_grib, dump_file);
660 if (header == dwd) wrtieee_header((int) len_grib, dump_file);
661 n_dump++;
662 }
663
664 if ((mode != INVENTORY && output_type != GRIB) || verbose > 1) {
665 /* decode numeric data */
666
667 if ((array = (float *) malloc(sizeof(float) * nxny)) == NULL) {
668 fprintf(stderr,"memory problems\n");
669 exit(8);
670 }
671
672 temp = int_power(10.0, - PDS_DecimalScale(pds));
673
674 BDS_unpack(array, bds, BMS_bitmap(bms), BDS_NumBits(bds), nxny,
675 temp*BDS_RefValue(bds),temp*int_power(2.0, BDS_BinScale(bds)));
676
677 if (verbose > 1) {
678 rmin = FLT_MAX;
679 rmax = -FLT_MAX;
680 for (i = 0; i < nxny; i++) {
681 if (fabs(array[i]-UNDEFINED) > 0.0001*UNDEFINED) {
682 rmin = min(rmin,array[i]);
683 rmax = max(rmax,array[i]);
684 }
685 }
686 printf(" min/max data %g %g num bits %d "
687 " BDS_Ref %g DecScale %d BinScale %d\n",
688 rmin, rmax, BDS_NumBits(bds), BDS_RefValue(bds),
689 PDS_DecimalScale(pds), BDS_BinScale(bds));
690 }
691
692 if (mode != INVENTORY && output_type != GRIB) {
693 /* dump code */
694 if (output_PDS_GDS == 1) {
695 /* insert code here */
696 if (output_type == BINARY || output_type == IEEE) {
697 /* write PDS */
698 i = PDS_LEN(pds) + 4;
699 if (header == simple && output_type == BINARY)
700 fwrite((void *) &i, sizeof(int), 1, dump_file);
701 if (header == simple && output_type == IEEE) wrtieee_header(i, dump_file);
702 fwrite((void *) "PDS ", 1, 4, dump_file);
703 fwrite((void *) pds, 1, i - 4, dump_file);
704 if (header == simple && output_type == BINARY)
705 fwrite((void *) &i, sizeof(int), 1, dump_file);
706 if (header == simple && output_type == IEEE) wrtieee_header(i, dump_file);
707
708 /* write GDS */
709 i = (gds) ? GDS_LEN(gds) + 4 : 4;
710 if (header == simple && output_type == BINARY)
711 fwrite((void *) &i, sizeof(int), 1, dump_file);
712 if (header == simple && output_type == IEEE) wrtieee_header(i, dump_file);
713 fwrite((void *) "GDS ", 1, 4, dump_file);
714 if (gds) fwrite((void *) gds, 1, i - 4, dump_file);
715 if (header == simple && output_type == BINARY)
716 fwrite((void *) &i, sizeof(int), 1, dump_file);
717 if (header == simple && output_type == IEEE) wrtieee_header(i, dump_file);
718 }
719 }
720
721 if (output_type == BINARY) {
722 i = nxny * sizeof(float);
723 if (header == simple) fwrite((void *) &i, sizeof(int), 1, dump_file);
724 fwrite((void *) array, sizeof(float), nxny, dump_file);
725 if (header == simple) fwrite((void *) &i, sizeof(int), 1, dump_file);
726 }
727 else if (output_type == IEEE) {
728 wrtieee(array, nxny, header, dump_file);
729 }
730 else if (output_type == TEXT) {
731 /* number of points in grid */
732 if (header == simple) {
733 if (nx <= 0 || ny <= 0 || nxny != nx*ny) {
734 fprintf(dump_file, "%ld %d\n", nxny, 1);
735 }
736 else {
737 fprintf(dump_file, "%d %d\n", nx, ny);
738 }
739 }
740 for (i = 0; i < nxny; i++) {
741 fprintf(dump_file,"%g\n", array[i]);
742 }
743 }
744 n_dump++;
745 }
746 free(array);
747 if (verbose > 0) printf("\n");
748 }
749
750 pos += len_grib;
751 count++;
752 }
753
754 if (mode != INVENTORY) {
755 if (header == dwd && output_type == GRIB) wrtieee_header(0, dump_file);
756 if (ferror(dump_file)) {
757 fprintf(stderr,"error writing %s\n",dump_file_name);
758 exit(8);
759 }
760 }
761 fclose(input);
762 return (return_code);
763 }
764
765 void print_pds(unsigned char *pds, int print_PDS, int print_PDS10, int verbose) {
766 int i, j;
767
768 j = PDS_LEN(pds);
769 if (verbose < 2) {
770 if (print_PDS && verbose < 2) {
771 printf(":PDS=");
772 for (i = 0; i < j; i++) {
773 printf("%2.2x", (int) pds[i]);
774 }
775 }
776 if (print_PDS10 && verbose < 2) {
777 printf(":PDS10=");
778 for (i = 0; i < j; i++) {
779 printf(" %d", (int) pds[i]);
780 }
781 }
782 }
783 else {
784 if (print_PDS) {
785 printf(" PDS(1..%d)=",j);
786 for (i = 0; i < j; i++) {
787 if (i % 20 == 0) printf("\n %4d:",i+1);
788 printf(" %3.2x", (int) pds[i]);
789 }
790 printf("\n");
791 }
792 if (print_PDS10) {
793 printf(" PDS10(1..%d)=",j);
794 for (i = 0; i < j; i++) {
795 if (i % 20 == 0) printf("\n %4d:",i+1);
796 printf(" %3d", (int) pds[i]);
797 }
798 printf("\n");
799 }
800 }
801 }
802
803 void print_gds(unsigned char *gds, int print_GDS, int print_GDS10, int verbose) {
804 int i, j;
805
806 j = GDS_LEN(gds);
807 if (verbose < 2) {
808 if (print_GDS && verbose < 2) {
809 printf(":GDS=");
810 for (i = 0; i < j; i++) {
811 printf("%2.2x", (int) gds[i]);
812 }
813 }
814 if (print_GDS10 && verbose < 2) {
815 printf(":GDS10=");
816 for (i = 0; i < j; i++) {
817 printf(" %d", (int) gds[i]);
818 }
819 }
820 }
821 else {
822 if (print_GDS) {
823 printf(" GDS(1..%d)=",j);
824 for (i = 0; i < j; i++) {
825 if (i % 20 == 0) printf("\n %4d:",i+1);
826 printf(" %3.2x", (int) gds[i]);
827 }
828 printf("\n");
829 }
830 if (print_GDS10) {
831 printf(" GDS10(1..%d)=",j);
832 for (i = 0; i < j; i++) {
833 if (i % 20 == 0) printf("\n %4d:",i+1);
834 printf(" %3d", (int) gds[i]);
835 }
836 printf("\n");
837 }
838 }
839 }
Mirror of the official repository at https://github.com/wrf-model/WRF