| blob | d8978b2619ea5a735bc1ebe660f023beed4ae79c |
1 subroutine compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
2 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
3 ! . . . .
4 ! SUBPROGRAM: compack
5 ! PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-06-21
6 !
7 ! ABSTRACT: This subroutine packs up a data field using a complex
8 ! packing algorithm as defined in the GRIB2 documention. It
9 ! supports GRIB2 complex packing templates with or without
10 ! spatial differences (i.e. DRTs 5.2 and 5.3).
11 ! It also fills in GRIB2 Data Representation Template 5.2 or 5.3
12 ! with the appropriate values.
13 !
14 ! PROGRAM HISTORY LOG:
15 ! 2000-06-21 Gilbert
16 !
17 ! USAGE: CALL compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
18 ! INPUT ARGUMENT LIST:
19 ! fld() - Contains the data values to pack
20 ! ndpts - The number of data values in array fld()
21 ! idrsnum - Data Representation Template number 5.N
22 ! Must equal 2 or 3.
23 ! idrstmpl - Contains the array of values for Data Representation
24 ! Template 5.2 or 5.3
25 ! (1) = Reference value - ignored on input
26 ! (2) = Binary Scale Factor
27 ! (3) = Decimal Scale Factor
28 ! .
29 ! .
30 ! (7) = Missing value management
31 ! (8) = Primary missing value
32 ! (9) = Secondary missing value
33 ! .
34 ! .
35 ! (17) = Order of Spatial Differencing ( 1 or 2 )
36 ! .
37 ! .
38 !
39 ! OUTPUT ARGUMENT LIST:
40 ! idrstmpl - Contains the array of values for Data Representation
41 ! Template 5.3
42 ! (1) = Reference value - set by compack routine.
43 ! (2) = Binary Scale Factor - unchanged from input
44 ! (3) = Decimal Scale Factor - unchanged from input
45 ! .
46 ! .
47 ! cpack - The packed data field (character*1 array)
48 ! lcpack - length of packed field cpack().
49 !
50 ! REMARKS: None
51 !
52 ! ATTRIBUTES:
53 ! LANGUAGE: XL Fortran 90
54 ! MACHINE: IBM SP
55 !
56 !$$$
58 integer,intent(in) :: ndpts,idrsnum
59 real,intent(in) :: fld(ndpts)
60 character(len=1),intent(out) :: cpack(*)
61 integer,intent(inout) :: idrstmpl(*)
62 integer,intent(out) :: lcpack
64 real(4) :: ref
65 integer(4) :: iref
66 integer,allocatable :: ifld(:)
67 integer,allocatable :: jmin(:),jmax(:),lbit(:)
68 integer,parameter :: zero=0
69 integer,allocatable :: gref(:),gwidth(:),glen(:)
70 integer :: glength,grpwidth
71 logical :: simple_alg = .false.
73 alog2=alog(2.0)
74 bscale=2.0**real(-idrstmpl(2))
75 dscale=10.0**real(idrstmpl(3))
76 !
77 ! Find max and min values in the data
78 !
79 rmax=fld(1)
80 rmin=fld(1)
81 do j=2,ndpts
82 if (fld(j).gt.rmax) rmax=fld(j)
83 if (fld(j).lt.rmin) rmin=fld(j)
84 enddo
85 !
86 ! If max and min values are not equal, pack up field.
87 ! If they are equal, we have a constant field, and the reference
88 ! value (rmin) is the value for each point in the field and
89 ! set nbits to 0.
90 !
91 if (rmin.ne.rmax) then
92 iofst=0
93 allocate(ifld(ndpts))
94 allocate(gref(ndpts))
95 allocate(gwidth(ndpts))
96 allocate(glen(ndpts))
97 !
98 ! Scale original data
99 !
100 if (idrstmpl(2).eq.0) then ! No binary scaling
101 imin=nint(rmin*dscale)
102 !imax=nint(rmax*dscale)
103 rmin=real(imin)
104 do j=1,ndpts
105 ifld(j)=nint(fld(j)*dscale)-imin
106 enddo
107 else ! Use binary scaling factor
108 rmin=rmin*dscale
109 !rmax=rmax*dscale
110 do j=1,ndpts
111 ifld(j)=nint(((fld(j)*dscale)-rmin)*bscale)
112 enddo
113 endif
114 !
115 ! Calculate Spatial differences, if using DRS Template 5.3
116 !
117 if (idrsnum.eq.3) then ! spatial differences
118 if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2
119 if (idrstmpl(17).eq.1) then ! first order
120 ival1=ifld(1)
121 do j=ndpts,2,-1
122 ifld(j)=ifld(j)-ifld(j-1)
123 enddo
124 ifld(1)=0
125 elseif (idrstmpl(17).eq.2) then ! second order
126 ival1=ifld(1)
127 ival2=ifld(2)
128 do j=ndpts,3,-1
129 ifld(j)=ifld(j)-(2*ifld(j-1))+ifld(j-2)
130 enddo
131 ifld(1)=0
132 ifld(2)=0
133 endif
134 !
135 ! subtract min value from spatial diff field
136 !
137 isd=idrstmpl(17)+1
138 minsd=minval(ifld(isd:ndpts))
139 do j=isd,ndpts
140 ifld(j)=ifld(j)-minsd
141 enddo
142 !
143 ! find num of bits need to store minsd and add 1 extra bit
144 ! to indicate sign
145 !
146 temp=alog(real(abs(minsd)+1))/alog2
147 nbitsd=ceiling(temp)+1
148 !
149 ! find num of bits need to store ifld(1) ( and ifld(2)
150 ! if using 2nd order differencing )
151 !
152 maxorig=ival1
153 if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2
154 temp=alog(real(maxorig+1))/alog2
155 nbitorig=ceiling(temp)+1
156 if (nbitorig.gt.nbitsd) nbitsd=nbitorig
157 ! increase number of bits to even multiple of 8 ( octet )
158 if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8))
159 !
160 ! Store extra spatial differencing info into the packed
161 ! data section.
162 !
163 if (nbitsd.ne.0) then
164 ! pack first original value
165 if (ival1.ge.0) then
166 call g2lib_sbyte(cpack,ival1,iofst,nbitsd)
167 iofst=iofst+nbitsd
168 else
169 call g2lib_sbyte(cpack,1,iofst,1)
170 iofst=iofst+1
171 call g2lib_sbyte(cpack,iabs(ival1),iofst,nbitsd-1)
172 iofst=iofst+nbitsd-1
173 endif
174 if (idrstmpl(17).eq.2) then
175 ! pack second original value
176 if (ival2.ge.0) then
177 call g2lib_sbyte(cpack,ival2,iofst,nbitsd)
178 iofst=iofst+nbitsd
179 else
180 call g2lib_sbyte(cpack,1,iofst,1)
181 iofst=iofst+1
182 call g2lib_sbyte(cpack,iabs(ival2),iofst,nbitsd-1)
183 iofst=iofst+nbitsd-1
184 endif
185 endif
186 ! pack overall min of spatial differences
187 if (minsd.ge.0) then
188 call g2lib_sbyte(cpack,minsd,iofst,nbitsd)
189 iofst=iofst+nbitsd
190 else
191 call g2lib_sbyte(cpack,1,iofst,1)
192 iofst=iofst+1
193 call g2lib_sbyte(cpack,iabs(minsd),iofst,nbitsd-1)
194 iofst=iofst+nbitsd-1
195 endif
196 endif
197 !print *,'SDp ',ival1,ival2,minsd,nbitsd
198 endif ! end of spatial diff section
199 !
200 ! Determine Groups to be used.
201 !
202 if ( simple_alg ) then
203 ! set group length to 10 : calculate number of groups
204 ! and length of last group
205 ngroups=ndpts/10
206 glen(1:ngroups)=10
207 itemp=mod(ndpts,10)
208 if (itemp.ne.0) then
209 ngroups=ngroups+1
210 glen(ngroups)=itemp
211 endif
212 else
213 ! Use Dr. Glahn's algorithm for determining grouping.
214 !
215 kfildo=6
216 minpk=10
217 inc=1
218 maxgrps=(ndpts/minpk)+1
219 allocate(jmin(maxgrps))
220 allocate(jmax(maxgrps))
221 allocate(lbit(maxgrps))
222 missopt=0
223 call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2,
224 & jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit,
225 & kbit,novref,lbitref,ier)
226 !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
227 do ng=1,ngroups
228 glen(ng)=glen(ng)+novref
229 enddo
230 deallocate(jmin)
231 deallocate(jmax)
232 deallocate(lbit)
233 endif
234 !
235 ! For each group, find the group's reference value
236 ! and the number of bits needed to hold the remaining values
237 !
238 n=1
239 do ng=1,ngroups
240 ! find max and min values of group
241 gref(ng)=ifld(n)
242 imax=ifld(n)
243 j=n+1
244 do lg=2,glen(ng)
245 if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j)
246 if (ifld(j).gt.imax) imax=ifld(j)
247 j=j+1
248 enddo
249 ! calc num of bits needed to hold data
250 if ( gref(ng).ne.imax ) then
251 temp=alog(real(imax-gref(ng)+1))/alog2
252 gwidth(ng)=ceiling(temp)
253 else
254 gwidth(ng)=0
255 endif
256 ! Subtract min from data
257 j=n
258 do lg=1,glen(ng)
259 ifld(j)=ifld(j)-gref(ng)
260 j=j+1
261 enddo
262 ! increment fld array counter
263 n=n+glen(ng)
264 enddo
265 !
266 ! Find max of the group references and calc num of bits needed
267 ! to pack each groups reference value, then
268 ! pack up group reference values
269 !
270 !write(77,*)'GREFS: ',(gref(j),j=1,ngroups)
271 igmax=maxval(gref(1:ngroups))
272 if (igmax.ne.0) then
273 temp=alog(real(igmax+1))/alog2
274 nbitsgref=ceiling(temp)
275 call g2lib_sbytes(cpack,gref,iofst,nbitsgref,0,ngroups)
276 itemp=nbitsgref*ngroups
277 iofst=iofst+itemp
278 ! Pad last octet with Zeros, if necessary,
279 if (mod(itemp,8).ne.0) then
280 left=8-mod(itemp,8)
281 call g2lib_sbyte(cpack,zero,iofst,left)
282 iofst=iofst+left
283 endif
284 else
285 nbitsgref=0
286 endif
287 !
288 ! Find max/min of the group widths and calc num of bits needed
289 ! to pack each groups width value, then
290 ! pack up group width values
291 !
292 !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
293 iwmax=maxval(gwidth(1:ngroups))
294 ngwidthref=minval(gwidth(1:ngroups))
295 if (iwmax.ne.ngwidthref) then
296 temp=alog(real(iwmax-ngwidthref+1))/alog2
297 nbitsgwidth=ceiling(temp)
298 do i=1,ngroups
299 gwidth(i)=gwidth(i)-ngwidthref
300 enddo
301 call g2lib_sbytes(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
302 itemp=nbitsgwidth*ngroups
303 iofst=iofst+itemp
304 ! Pad last octet with Zeros, if necessary,
305 if (mod(itemp,8).ne.0) then
306 left=8-mod(itemp,8)
307 call g2lib_sbyte(cpack,zero,iofst,left)
308 iofst=iofst+left
309 endif
310 else
311 nbitsgwidth=0
312 gwidth(1:ngroups)=0
313 endif
314 !
315 ! Find max/min of the group lengths and calc num of bits needed
316 ! to pack each groups length value, then
317 ! pack up group length values
318 !
319 !write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
320 ilmax=maxval(glen(1:ngroups-1))
321 nglenref=minval(glen(1:ngroups-1))
322 nglenlast=glen(ngroups)
323 if (ilmax.ne.nglenref) then
324 temp=alog(real(ilmax-nglenref+1))/alog2
325 nbitsglen=ceiling(temp)
326 do i=1,ngroups-1
327 glen(i)=glen(i)-nglenref
328 enddo
329 call g2lib_sbytes(cpack,glen,iofst,nbitsglen,0,ngroups)
330 itemp=nbitsglen*ngroups
331 iofst=iofst+itemp
332 ! Pad last octet with Zeros, if necessary,
333 if (mod(itemp,8).ne.0) then
334 left=8-mod(itemp,8)
335 call g2lib_sbyte(cpack,zero,iofst,left)
336 iofst=iofst+left
337 endif
338 else
339 nbitsglen=0
340 glen(1:ngroups)=0
341 endif
342 !
343 ! For each group, pack data values
344 !
345 !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
346 n=1
347 ij=0
348 do ng=1,ngroups
349 glength=glen(ng)+nglenref
350 if (ng.eq.ngroups ) glength=nglenlast
351 grpwidth=gwidth(ng)+ngwidthref
352 !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
353 if ( grpwidth.ne.0 ) then
354 call g2lib_sbytes(cpack,ifld(n),iofst,grpwidth,0,glength)
355 iofst=iofst+(grpwidth*glength)
356 endif
357 do kk=1,glength
358 ij=ij+1
359 !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
360 enddo
361 n=n+glength
362 enddo
363 ! Pad last octet with Zeros, if necessary,
364 if (mod(iofst,8).ne.0) then
365 left=8-mod(iofst,8)
366 call g2lib_sbyte(cpack,zero,iofst,left)
367 iofst=iofst+left
368 endif
369 lcpack=iofst/8
370 !
371 if ( allocated(ifld) ) deallocate(ifld)
372 if ( allocated(gref) ) deallocate(gref)
373 if ( allocated(gwidth) ) deallocate(gwidth)
374 if ( allocated(glen) ) deallocate(glen)
375 else ! Constant field ( max = min )
376 nbits=0
377 lcpack=0
378 nbitsgref=0
379 ngroups=0
380 endif
382 !
383 ! Fill in ref value and number of bits in Template 5.2
384 !
385 call mkieee(rmin,ref,1) ! ensure reference value is IEEE format
386 ! call g2lib_gbyte(ref,idrstmpl(1),0,32)
387 iref=transfer(ref,iref)
388 idrstmpl(1)=iref
389 idrstmpl(4)=nbitsgref
390 idrstmpl(5)=0 ! original data were reals
391 idrstmpl(6)=1 ! general group splitting
392 idrstmpl(7)=0 ! No internal missing values
393 idrstmpl(8)=0 ! Primary missing value
394 idrstmpl(9)=0 ! secondary missing value
395 idrstmpl(10)=ngroups ! Number of groups
396 idrstmpl(11)=ngwidthref ! reference for group widths
397 idrstmpl(12)=nbitsgwidth ! num bits used for group widths
398 idrstmpl(13)=nglenref ! Reference for group lengths
399 idrstmpl(14)=1 ! length increment for group lengths
400 idrstmpl(15)=nglenlast ! True length of last group
401 idrstmpl(16)=nbitsglen ! num bits used for group lengths
402 if (idrsnum.eq.3) then
403 idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial
404 ! differencing values
405 endif
407 return
408 end