source: trunk/WRF.COMMON/WRFV3/external/io_grib1/WGRIB/wgrib_main.c @ 2759

Last change on this file since 2759 was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 29.3 KB
Line 
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
104enum Def_NCEP_Table def_ncep_table = DEF_T62_NCEP_TABLE;
105int minute = 0;
106int ncep_ens = 0;
107
108int 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
765void 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
803void 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}
Note: See TracBrowser for help on using the repository browser.