-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdpugrb.f
95 lines (95 loc) · 2.27 KB
/
dpugrb.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
SUBROUTINE DP_UGRB ( idata, kxky, nbits, qmin, scale, misflg,
+ grid, iret )
C************************************************************************
C* DP_UGRB *
C* *
C* This subroutine unpacks a grid into the GEMPAK GRIB format. *
C* The packing and unpacking equations are: *
C* *
C* IDATA = NINT ( ( GRID - QMIN ) / SCALE ) *
C* GRID = QMIN + IDATA * SCALE *
C* *
C* DP_UGRB ( IDATA, KXKY, NBITS, QMIN, SCALE, MISFLG, GRID, IRET ) *
C* *
C* Input parameters: *
C* IDATA (*) INTEGER Packed data *
C* KXKY INTEGER Number of grid points *
C* NBITS INTEGER Number of bits *
C* QMIN REAL Minimum value of grid *
C* SCALE REAL Scaling factor *
C* MISFLG LOGICAL Missing data flag *
C* *
C* Output parameters: *
C* GRID (KXKY) REAL Grid data *
C* IRET INTEGER Return code *
C* 0 = normal return *
C* -10 = NBITS invalid *
C* -12 = invalid scale *
C** *
C* Log: *
C* M. desJardins/GSFC 3/89 *
C************************************************************************
INCLUDE 'GEMPRM.PRM'
C*
REAL grid (*)
INTEGER idata (*)
LOGICAL misflg
C------------------------------------------------------------------------
iret = 0
C
C* Check for valid input.
C
IF ( ( nbits .le. 1 ) .or. ( nbits .gt. 31 ) ) THEN
iret = -10
RETURN
END IF
IF ( scale .eq. 0. ) THEN
iret = -12
RETURN
END IF
C
C* Compute missing data value.
C
imax = 2 ** nbits - 1
C
C* Retrieve data points from buffer.
C
iword = 1
ibit = 1
DO i = 1, kxky
C
C* Get the integer from the buffer.
C
jshft = nbits + ibit - 33
idat = 0
idat = ISHFT ( idata (iword), jshft )
idat = IAND ( idat, imax )
C
C* Check to see if packed integer overflows into next word.
C
IF ( jshft .gt. 0 ) THEN
jshft = jshft - 32
idat2 = 0
idat2 = ISHFT ( idata (iword+1), jshft )
idat = IOR ( idat, idat2 )
END IF
C
C* Compute value of word.
C
IF ( ( idat .eq. imax ) .and. misflg ) THEN
grid (i) = RMISSD
ELSE
grid (i) = qmin + idat * scale
END IF
C
C* Set location for next word.
C
ibit = ibit + nbits
IF ( ibit .gt. 32 ) THEN
ibit = ibit - 32
iword = iword + 1
END IF
END DO
C*
RETURN
END