1 #include "grib_int.h"
2
3
confp3(double pval,int * kexp,int * kmant,int kbits,int kround)4 void confp3(double pval, int *kexp, int *kmant, int kbits, int kround)
5 {
6 /*
7
8 Purpose:
9 --------
10
11 Convert floating point number from machine
12 representation to GRIB representation.
13
14 Input Parameters:
15 -----------------
16
17 pval - Floating point number to be converted.
18 kbits - Number of bits in computer word.
19 kround - Conversion type.
20 0 , Closest number in GRIB format less than
21 original number.
22 1 , Closest number in GRIB format to the
23 original number (equal to, greater than or
24 less than original number).
25
26 Output Parameters:
27 ------------------
28
29 kexp - 8 Bit signed exponent.
30 kmant - 24 Bit mantissa.
31
32 Method:
33 -------
34
35 Floating point number represented as 8 bit signed
36 exponent and 24 bit mantissa in integer values.
37
38 Externals.
39 ----------
40
41 decfp2 - Decode from IBM floating point format.
42
43 Reference:
44 ----------
45
46 WMO Manual on Codes re GRIB representation.
47
48 Comments:
49 ---------
50
51 Routine aborts if an invalid conversion type parameter
52 is used or if a 24 bit mantissa is not produced.
53
54 Author:
55 -------
56
57 John Hennessy ECMWF 18.06.91
58
59 Modifications:
60 --------------
61
62 Uwe Schulzweida MPIfM 01/04/2001
63
64 Convert to C from EMOS library version 130
65
66 Uwe Schulzweida MPIfM 02/08/2002
67
68 - speed up by factor 1.6 on NEC SX6
69 - replace 1.0 / pow(16.0, (double)(iexp - 70)) by rpow16m70tab[iexp]
70 */
71
72 // extern int CGRIBEX_Debug;
73
74 /* ----------------------------------------------------------------- */
75 /* Section 1 . Initialise */
76 /* ----------------------------------------------------------------- */
77
78 // Check conversion type parameter.
79
80 int iround = kround;
81 if ( iround != 0 && iround != 1 )
82 {
83 Error("Invalid conversion type = %d", iround);
84
85 // If not aborting, arbitrarily set rounding to 'up'.
86 iround = 1;
87 }
88
89 /* ----------------------------------------------------------------- */
90 /* Section 2 . Convert value of zero. */
91 /* ----------------------------------------------------------------- */
92
93 if (fabs(pval) <= 0)
94 {
95 *kexp = 0;
96 *kmant = 0;
97 goto LABEL900;
98 }
99
100 /* ----------------------------------------------------------------- */
101 /* Section 3 . Convert other values. */
102 /* ----------------------------------------------------------------- */
103 {
104 const double zeps = (kbits != 32) ? 1.0e-12 : 1.0e-8;
105 double zref = pval;
106
107 // Sign of value.
108 const int isign = (zref >= 0.0) ? 0 : 128;
109 zref = fabs(zref);
110
111 // Exponent.
112 int iexp = (int) (log(zref)/log(16.0) + 65.0 + zeps);
113
114 // only ANSI C99 has log2
115 // iexp = (int) (log2(zref) * 0.25 + 65.0 + zeps);
116
117 if ( iexp < 0 ) iexp = 0;
118 if ( iexp > 127 ) iexp = 127;
119
120 // double rpowref = zref / pow(16.0, (double)(iexp - 70));
121 double rpowref = ldexp(zref, 4 * -(iexp - 70));
122
123 // Mantissa.
124 if ( iround == 0 )
125 {
126 /* Closest number in GRIB format less than original number. */
127 /* Truncate for positive numbers. */
128 /* Round up for negative numbers. */
129 *kmant = (isign == 0) ? (int)rpowref : (int)lround(rpowref + 0.5);
130 }
131 else
132 {
133 /* Closest number in GRIB format to the original number */
134 /* (equal to, greater than or less than original number). */
135 *kmant = (int)lround(rpowref);
136 }
137
138 /* Check that mantissa value does not exceed 24 bits. */
139 /* If it does, adjust the exponent upwards and recalculate the mantissa. */
140 /* 16777215 = 2**24 - 1 */
141 if ( *kmant > 16777215 )
142 {
143
144 LABEL350:
145
146 ++iexp;
147
148 // Check for exponent overflow during adjustment
149 if ( iexp > 127 )
150 {
151 Message("Exponent overflow");
152 Message("Original number = %30.20f", pval);
153 Message("Sign = %3d, Exponent = %3d, Mantissa = %12d", isign, iexp, *kmant);
154
155 Error("Exponent overflow");
156
157 // If not aborting, arbitrarily set value to zero
158 Message("Value arbitrarily set to zero.");
159 *kexp = 0;
160 *kmant = 0;
161 goto LABEL900;
162 }
163
164 rpowref = ldexp(zref, 4 * -(iexp - 70));
165
166 if ( iround == 0 )
167 {
168 /* Closest number in GRIB format less than original number. */
169 /* Truncate for positive numbers. */
170 /* Round up for negative numbers. */
171 *kmant = (isign == 0) ? (int)rpowref : (int)lround(rpowref + 0.5);
172 }
173 else
174 {
175 /* Closest number in GRIB format to the original number */
176 /* (equal to, greater or less than original number). */
177 *kmant = (int)lround(rpowref);
178 }
179
180 // Repeat calculation (with modified exponent) if still have mantissa overflow.
181 if ( *kmant > 16777215 ) goto LABEL350;
182 }
183
184 // Add sign bit to exponent.
185 *kexp = iexp + isign;
186 }
187
188 /* ----------------------------------------------------------------- */
189 /* Section 9. Return */
190 /* ----------------------------------------------------------------- */
191
192 LABEL900:
193 /*
194 if ( CGRIBEX_Debug )
195 {
196 double zval;
197
198 Message("Conversion type parameter = %4d", kround);
199 Message("Original number = %30.20f", pval);
200
201 zval = decfp2(*kexp, *kmant);
202
203 Message("Converted to %30.20f", zval);
204 Message("Sign = %3d, Exponent = %3d, Mantissa = %12d", isign, iexp, *kmant);
205 }
206 */
207 return;
208 } /* confp3 */
209