source: lmdz_wrf/trunk/WRFV3/external/io_grib1/WGRIB/BDSunpk.c @ 14

Last change on this file since 14 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 3.4 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <stddef.h>
4#include <limits.h>
5#include "grib.h"
6#include "pds4.h"
7#include "bms.h"
8#include "bds.h"
9
10/* 1996                         wesley ebisuzaki
11 *
12 * Unpack BDS section
13 *
14 * input: *bits, pointer to packed integer data
15 *        *bitmap, pointer to bitmap (undefined data), NULL if none
16 *        n_bits, number of bits per packed integer
17 *        n, number of data points (includes undefined data)
18 *        ref, scale: flt[] = ref + scale*packed_int
19 * output: *flt, pointer to output array
20 *        undefined values filled with UNDEFINED
21 *
22 * note: code assumes an integer > 32 bits
23 *
24 * 7/98 v1.2.1 fix bug for bitmaps and nbit >= 25 found by Larry Brasfield
25 * 2/01 v1.2.2 changed jj from long int to double
26 * 3/02 v1.2.3 added unpacking extensions for spectral data
27 *             Luis Kornblueh, MPIfM
28 */
29
30static unsigned int mask[] = {0,1,3,7,15,31,63,127,255};
31static unsigned int map_masks[8] = {128, 64, 32, 16, 8, 4, 2, 1};
32static double shift[9] = {1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0};
33
34void BDS_unpack(float *flt, unsigned char *bds, unsigned char *bitmap,
35        int n_bits, int n, double ref, double scale) {
36
37    unsigned char *bits;
38
39    int i, mask_idx, t_bits, c_bits, j_bits;
40    unsigned int j, map_mask, tbits, jmask, bbits;
41    double jj;
42
43    if (BDS_Harmonic(bds)) {
44        bits = bds + 15;
45        /* fill in global mean */
46        *flt++ = BDS_Harmonic_RefValue(bds);
47        n -= 1; 
48    }
49    else {
50        bits = bds + 11; 
51    }
52
53    tbits = bbits = 0;
54
55    /* assume integer has 32+ bits */
56    if (n_bits <= 25) {
57        jmask = (1 << n_bits) - 1;
58        t_bits = 0;
59
60        if (bitmap) {
61            for (i = 0; i < n; i++) {
62                /* check bitmap */
63                mask_idx = i & 7;
64                if (mask_idx == 0) bbits = *bitmap++;
65                if ((bbits & map_masks[mask_idx]) == 0) {
66                    *flt++ = UNDEFINED;
67                    continue;
68                }
69
70                while (t_bits < n_bits) {
71                    tbits = (tbits * 256) + *bits++;
72                    t_bits += 8;
73                }
74                t_bits -= n_bits;
75                j = (tbits >> t_bits) & jmask;
76                *flt++ = ref + scale*j;
77            }
78        }
79        else {
80            for (i = 0; i < n; i++) {
81                while (t_bits < n_bits) {
82                    tbits = (tbits * 256) + *bits++;
83                    t_bits += 8;
84                }
85                t_bits -= n_bits;
86                flt[i] = (tbits >> t_bits) & jmask;
87            }
88            /* at least this vectorizes :) */
89            for (i = 0; i < n; i++) {
90                flt[i] = ref + scale*flt[i];
91            }
92        }
93    }
94    else {
95        /* older unoptimized code, not often used */
96        c_bits = 8;
97        map_mask = 128;
98        while (n-- > 0) {
99            if (bitmap) {
100                j = (*bitmap & map_mask);
101                if ((map_mask >>= 1) == 0) {
102                    map_mask = 128;
103                    bitmap++;
104                }
105                if (j == 0) {
106                    *flt++ = UNDEFINED;
107                    continue;
108                }
109            }
110
111            jj = 0.0;
112            j_bits = n_bits;
113            while (c_bits <= j_bits) {
114                if (c_bits == 8) {
115                    jj = jj * 256.0  + (double) (*bits++);
116                    j_bits -= 8;
117                }
118                else {
119                    jj = (jj * shift[c_bits]) + (double) (*bits & mask[c_bits]);
120                    bits++;
121                    j_bits -= c_bits;
122                    c_bits = 8;
123                }
124            }
125            if (j_bits) {
126                c_bits -= j_bits;
127                jj = (jj * shift[j_bits]) + (double) ((*bits >> c_bits) & mask[j_bits]);
128            }
129            *flt++ = ref + scale*jj;
130        }
131    }
132    return;
133}
Note: See TracBrowser for help on using the repository browser.