1 #include <stdio.h>
2 #include <string.h>
3 #include <stdlib.h>
4 #include <stdbool.h>
5 #include <unistd.h>
6 #include <math.h>
7
8 int timer_new(const char *text);
9 void timer_report(void);
10 void timer_start(int it);
11 void timer_stop(int it);
12
13 //#undef _GET_IBM_COUNTER
14
15 #ifdef _GET_IBM_COUNTER
16 #include <libhpc.h>
17 #endif
18
19 #include "cgribex.h"
20
21 #define NINT(x) ((x) < 0 ? (int)((x)-.5) : (int)((x)+.5))
22
23 #define NLEV 1 // Number of levels
24 #define WITH_IO 1
25
26 static
cgribexDefaultSec0(int * isec0)27 void cgribexDefaultSec0(int *isec0)
28 {
29 ISEC0_GRIB_Len = 0;
30 ISEC0_GRIB_Version = 0;
31 }
32
33 static
cgribexDefaultSec1(int * isec1)34 void cgribexDefaultSec1(int *isec1)
35 {
36 ISEC1_CenterID = 0;
37 ISEC1_SubCenterID = 0;
38 ISEC1_LocalFLag = 0;
39 }
40
41 static
cgribexDefaultSec4(int * isec4)42 void cgribexDefaultSec4(int *isec4)
43 {
44 for ( long i = 2; i <= 10; ++i ) isec4[i] = 0;
45 }
46
47 static
cgribexDefGridLonLat(int * isec1,int * isec2,int * isec4,int nlon,int nlat,double di)48 void cgribexDefGridLonLat(int *isec1, int *isec2, int *isec4, int nlon, int nlat, double di)
49 {
50 memset(isec2, 0, 16*sizeof(int));
51
52 double xfirst = 0;
53 double xlast = 360. - di/2;
54 double yfirst = -90. + di/2;
55 double ylast = 90. - di/2;
56
57 ISEC1_Sec2Or3Flag = 128;
58 ISEC1_GridDefinition = 255;
59
60 ISEC2_GridType = GRIB1_GTYPE_LATLON;
61 ISEC2_NumLon = nlon;
62 ISEC2_NumLat = nlat;
63 ISEC2_FirstLat = NINT(yfirst*1000);
64 ISEC2_LastLat = NINT(ylast*1000);
65 ISEC2_LatIncr = NINT(di*1000);
66 ISEC2_FirstLon = NINT(xfirst*1000);
67 ISEC2_LastLon = NINT(xlast*1000);
68 ISEC2_LonIncr = NINT(di*1000);
69
70 ISEC2_Reduced = 0;
71 ISEC2_ResFlag = 128;
72 ISEC2_ScanFlag = 128;
73 ISEC2_ScanFlag += 64;
74 }
75
76 static
cgribexDefLevelSurface(int * isec1,int * isec2,double * fsec2)77 void cgribexDefLevelSurface(int *isec1, int *isec2, double *fsec2)
78 {
79 ISEC2_NumVCP = 0;
80 ISEC1_LevelType = GRIB1_LTYPE_SURFACE;
81 ISEC1_Level1 = 0;
82 ISEC1_Level2 = 0;
83 }
84
85 static
cgribexDefTime(int * isec1,int time)86 void cgribexDefTime(int *isec1, int time)
87 {
88 int year = 2000, month = 1, day = 1;
89
90 int hour = time/10000;
91 int minute = (time - hour*10000) / 100;
92
93 int century = year / 100;
94
95 ISEC1_Year = year - century*100;
96
97 if ( year < 0 )
98 {
99 century = -century;
100 ISEC1_Year = -ISEC1_Year;
101 }
102
103 if ( ISEC1_Year == 0 )
104 {
105 century -= 1;
106 ISEC1_Year = 100;
107 }
108
109 century += 1;
110 if ( year < 0 ) century = -century;
111
112 ISEC1_Month = month;
113 ISEC1_Day = day;
114 ISEC1_Hour = hour;
115 ISEC1_Minute = minute;
116
117 ISEC1_Century = century;
118
119 ISEC1_TimeUnit = ISEC1_TABLE4_MINUTE;
120 ISEC1_TimeRange = 10;
121
122 ISEC1_TimePeriod1 = 0;
123 ISEC1_TimePeriod2 = 0;
124 ISEC1_AvgNum = 0;
125 ISEC1_AvgMiss = 0;
126 ISEC1_DecScaleFactor = 0;
127 }
128
usage(void)129 void usage(void)
130 {
131 fprintf(stderr,"\ngrib_encode [-b <bits per sample> -d <grid direction increment> -n <num recs> -z -t] <filename>\n\n");
132 exit(-1);
133 }
134
main(int argc,char * argv[])135 int main(int argc, char *argv[])
136 {
137 int isec0[2], isec1[4096], isec2[4096], isec3[2], isec4[512];
138 double fsec2[512], fsec3[2];
139 int nbits = 16;
140 int c, iret, iword;
141 int fileID;
142 int nts = 100;
143 double di = 0.1; // grid direction increment
144 bool lzip = false;
145 bool gribTimer = false;
146
147 while ((c = getopt(argc, argv, "b:d:n:azt")) != -1)
148 {
149 switch (c)
150 {
151 case 'a':
152 case 'z':
153 lzip = true;
154 break;
155 case 'b':
156 nbits = atoi(optarg);
157 break;
158 case 'd':
159 di = atof(optarg);
160 break;
161 case 'n':
162 nts = atoi(optarg);
163 break;
164 case 't':
165 gribTimer = true;
166 break;
167 case '?':
168 default:
169 usage();
170 }
171 }
172 argc -= optind;
173 argv += optind;
174
175 if (nbits < 1) usage();
176 if ((argc == 0) | (argc > 1)) usage();
177
178 char *filename = (char *) malloc(sizeof(char)*(strlen(argv[0])+1));
179 strcpy(filename, argv[0]);
180
181 #ifdef _GET_IBM_COUNTER
182 hpmInit(0,"main");
183 #endif
184
185 int nlat = (int) (180./di);
186 int nlon = nlat*2;
187 int gridsize = nlon*nlat;
188
189 double *data = (double *) malloc(gridsize*sizeof(double));
190 for ( int i = 0; i < gridsize; ++i )
191 {
192 // data[i] = 273. + 50*(rand()/(RAND_MAX+1.0));
193 int ilat = (i%gridsize)/ nlon;
194 int ilon = i%nlon;
195 data[i] = cos(2.0 * M_PI * (0 + ilon*di)/360)
196 * sin(2.0 * M_PI * (-90 + di/2 + ilat*di)/180);
197 }
198
199 memset(isec1, 0, 32*sizeof(int));
200 fsec2[0] = 0; fsec2[1] = 0;
201
202 cgribexDefaultSec0(isec0);
203 cgribexDefaultSec1(isec1);
204 cgribexDefaultSec4(isec4);
205
206 ISEC1_CodeTable = 201;
207 ISEC1_Parameter = 1;
208
209 cgribexDefGridLonLat(isec1, isec2, isec4, nlon, nlat, di);
210 cgribexDefLevelSurface(isec1, isec2, fsec2);
211
212 if ( nbits > 0 ) isec4[1] = nbits;
213
214 ISEC4_NumValues = gridsize;
215 ISEC4_NumBits = nbits;
216
217 size_t gribbuffersize = isec4[0]*isec4[1]/8 + isec4[0]/8 + 3000; // datasize + mask + header
218 long gribsize = gribbuffersize / sizeof(int);
219 unsigned char *gribbuffer = (unsigned char *) malloc(gribbuffersize);
220 size_t buffersize = 0;
221 unsigned char *buffer = NULL;
222
223 int timer_total, timer_write, timer_encode, timer_szip;
224 if ( gribTimer )
225 {
226 timer_total = timer_new("total");
227 timer_write = timer_new("write");
228 timer_encode = timer_new("encode");
229 timer_szip = timer_new("szip");
230 }
231
232 if ( gribTimer ) timer_start(timer_total);
233
234 if ( WITH_IO ) fileID = gribOpen(filename, "w");
235
236 for ( int i = 0; i < nts; ++i )
237 {
238 cgribexDefTime(isec1, 120000+i*100);
239
240 if ( gribTimer ) timer_start(timer_encode);
241 gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, (double*) data,
242 gridsize, (int *) gribbuffer, gribsize, &iword, "C", &iret);
243 if ( gribTimer ) timer_stop(timer_encode);
244
245 if ( iret )
246 {
247 fprintf(stderr, "Problem during GRIB encode (errno = %d)!\n", iret);
248 break;
249 }
250
251 size_t nbytes = iword*sizeof(int);
252
253 if ( lzip )
254 {
255 if ( buffersize < gribbuffersize )
256 {
257 buffersize = gribbuffersize;
258 buffer = (unsigned char *) realloc(buffer, buffersize);
259 }
260
261 if ( gribTimer ) timer_start(timer_szip);
262 nbytes = gribZip(gribbuffer, gribbuffersize, buffer, buffersize);
263 if ( gribTimer ) timer_stop(timer_szip);
264 }
265
266 if ( gribTimer ) timer_start(timer_write);
267 if ( WITH_IO ) (void) gribWrite(fileID, gribbuffer, nbytes);
268 if ( gribTimer ) timer_stop(timer_write);
269 }
270
271 if ( WITH_IO ) gribClose(fileID);
272
273 free(gribbuffer);
274 free(data);
275 free(filename);
276
277 #ifdef _GET_IBM_COUNTER
278 hpmTerminate(0);
279 #endif
280
281 if ( gribTimer ) timer_stop(timer_total);
282 if ( gribTimer ) timer_report();
283
284 return 0;
285 }
286