| blob | 77816e35b5d2c6245859cdbe31b7d77b304c647c |
1 subroutine misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
2 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
3 ! . . . .
4 ! SUBPROGRAM: misspack
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 ! This version assumes that Missing Value Management is being used and that
14 ! 1 or 2 missing values appear in the data.
15 !
16 ! PROGRAM HISTORY LOG:
17 ! 2000-06-21 Gilbert
18 ! 2004-12-29 Gilbert - Corrected bug when encoding secondary missing values.
19 !
20 ! USAGE: CALL misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
21 ! INPUT ARGUMENT LIST:
22 ! fld() - Contains the data values to pack
23 ! ndpts - The number of data values in array fld()
24 ! idrsnum - Data Representation Template number 5.N
25 ! Must equal 2 or 3.
26 ! idrstmpl - Contains the array of values for Data Representation
27 ! Template 5.2 or 5.3
28 ! (1) = Reference value - ignored on input
29 ! (2) = Binary Scale Factor
30 ! (3) = Decimal Scale Factor
31 ! .
32 ! .
33 ! (7) = Missing value management
34 ! (8) = Primary missing value
35 ! (9) = Secondary missing value
36 ! .
37 ! .
38 ! (17) = Order of Spatial Differencing ( 1 or 2 )
39 ! .
40 ! .
41 !
42 ! OUTPUT ARGUMENT LIST:
43 ! idrstmpl - Contains the array of values for Data Representation
44 ! Template 5.3
45 ! (1) = Reference value - set by misspack routine.
46 ! (2) = Binary Scale Factor - unchanged from input
47 ! (3) = Decimal Scale Factor - unchanged from input
48 ! .
49 ! .
50 ! cpack - The packed data field (character*1 array)
51 ! lcpack - length of packed field cpack().
52 !
53 ! REMARKS: None
54 !
55 ! ATTRIBUTES:
56 ! LANGUAGE: XL Fortran 90
57 ! MACHINE: IBM SP
58 !
59 !$$$
61 integer,intent(in) :: ndpts,idrsnum
62 real,intent(in) :: fld(ndpts)
63 character(len=1),intent(out) :: cpack(*)
64 integer,intent(inout) :: idrstmpl(*)
65 integer,intent(out) :: lcpack
67 real(4) :: ref
68 integer(4) :: iref
69 integer,allocatable :: ifld(:),ifldmiss(:),jfld(:)
70 integer,allocatable :: jmin(:),jmax(:),lbit(:)
71 integer,parameter :: zero=0
72 integer,allocatable :: gref(:),gwidth(:),glen(:)
73 integer :: glength,grpwidth
74 logical :: simple_alg = .false.
76 alog2=alog(2.0)
77 bscale=2.0**real(-idrstmpl(2))
78 dscale=10.0**real(idrstmpl(3))
79 missopt=idrstmpl(7)
80 if ( missopt.ne.1 .AND. missopt.ne.2 ) then
81 print *,'misspack: Unrecognized option.'
82 lcpack=-1
83 return
84 else ! Get missing values
85 call rdieee(idrstmpl(8),rmissp,1)
86 if (missopt.eq.2) call rdieee(idrstmpl(9),rmisss,1)
87 endif
88 !
89 ! Find min value of non-missing values in the data,
90 ! AND set up missing value mapping of the field.
91 !
92 allocate(ifldmiss(ndpts))
93 rmin=huge(rmin)
94 if ( missopt .eq. 1 ) then ! Primary missing value only
95 do j=1,ndpts
96 if (fld(j).eq.rmissp) then
97 ifldmiss(j)=1
98 else
99 ifldmiss(j)=0
100 if (fld(j).lt.rmin) rmin=fld(j)
101 endif
102 enddo
103 endif
104 if ( missopt .eq. 2 ) then ! Primary and secondary missing values
105 do j=1,ndpts
106 if (fld(j).eq.rmissp) then
107 ifldmiss(j)=1
108 elseif (fld(j).eq.rmisss) then
109 ifldmiss(j)=2
110 else
111 ifldmiss(j)=0
112 if (fld(j).lt.rmin) rmin=fld(j)
113 endif
114 enddo
115 endif
116 !
117 ! Allocate work arrays:
118 ! Note: -ifldmiss(j),j=1,ndpts is a map of original field indicating
119 ! which of the original data values
120 ! are primary missing (1), sencondary missing (2) or non-missing (0).
121 ! -jfld(j),j=1,nonmiss is a subarray of just the non-missing values from
122 ! the original field.
123 !
124 !if (rmin.ne.rmax) then
125 iofst=0
126 allocate(ifld(ndpts))
127 allocate(jfld(ndpts))
128 allocate(gref(ndpts))
129 allocate(gwidth(ndpts))
130 allocate(glen(ndpts))
131 !
132 ! Scale original data
133 !
134 nonmiss=0
135 if (idrstmpl(2).eq.0) then ! No binary scaling
136 imin=nint(rmin*dscale)
137 !imax=nint(rmax*dscale)
138 rmin=real(imin)
139 do j=1,ndpts
140 if (ifldmiss(j).eq.0) then
141 nonmiss=nonmiss+1
142 jfld(nonmiss)=nint(fld(j)*dscale)-imin
143 endif
144 enddo
145 else ! Use binary scaling factor
146 rmin=rmin*dscale
147 !rmax=rmax*dscale
148 do j=1,ndpts
149 if (ifldmiss(j).eq.0) then
150 nonmiss=nonmiss+1
151 jfld(nonmiss)=nint(((fld(j)*dscale)-rmin)*bscale)
152 endif
153 enddo
154 endif
155 !
156 ! Calculate Spatial differences, if using DRS Template 5.3
157 !
158 if (idrsnum.eq.3) then ! spatial differences
159 if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2
160 if (idrstmpl(17).eq.1) then ! first order
161 ival1=jfld(1)
162 do j=nonmiss,2,-1
163 jfld(j)=jfld(j)-jfld(j-1)
164 enddo
165 jfld(1)=0
166 elseif (idrstmpl(17).eq.2) then ! second order
167 ival1=jfld(1)
168 ival2=jfld(2)
169 do j=nonmiss,3,-1
170 jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2)
171 enddo
172 jfld(1)=0
173 jfld(2)=0
174 endif
175 !
176 ! subtract min value from spatial diff field
177 !
178 isd=idrstmpl(17)+1
179 minsd=minval(jfld(isd:nonmiss))
180 do j=isd,nonmiss
181 jfld(j)=jfld(j)-minsd
182 enddo
183 !
184 ! find num of bits need to store minsd and add 1 extra bit
185 ! to indicate sign
186 !
187 temp=alog(real(abs(minsd)+1))/alog2
188 nbitsd=ceiling(temp)+1
189 !
190 ! find num of bits need to store ifld(1) ( and ifld(2)
191 ! if using 2nd order differencing )
192 !
193 maxorig=ival1
194 if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2
195 temp=alog(real(maxorig+1))/alog2
196 nbitorig=ceiling(temp)+1
197 if (nbitorig.gt.nbitsd) nbitsd=nbitorig
198 ! increase number of bits to even multiple of 8 ( octet )
199 if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8))
200 !
201 ! Store extra spatial differencing info into the packed
202 ! data section.
203 !
204 if (nbitsd.ne.0) then
205 ! pack first original value
206 if (ival1.ge.0) then
207 call g2lib_sbyte(cpack,ival1,iofst,nbitsd)
208 iofst=iofst+nbitsd
209 else
210 call g2lib_sbyte(cpack,1,iofst,1)
211 iofst=iofst+1
212 call g2lib_sbyte(cpack,iabs(ival1),iofst,nbitsd-1)
213 iofst=iofst+nbitsd-1
214 endif
215 if (idrstmpl(17).eq.2) then
216 ! pack second original value
217 if (ival2.ge.0) then
218 call g2lib_sbyte(cpack,ival2,iofst,nbitsd)
219 iofst=iofst+nbitsd
220 else
221 call g2lib_sbyte(cpack,1,iofst,1)
222 iofst=iofst+1
223 call g2lib_sbyte(cpack,iabs(ival2),iofst,nbitsd-1)
224 iofst=iofst+nbitsd-1
225 endif
226 endif
227 ! pack overall min of spatial differences
228 if (minsd.ge.0) then
229 call g2lib_sbyte(cpack,minsd,iofst,nbitsd)
230 iofst=iofst+nbitsd
231 else
232 call g2lib_sbyte(cpack,1,iofst,1)
233 iofst=iofst+1
234 call g2lib_sbyte(cpack,iabs(minsd),iofst,nbitsd-1)
235 iofst=iofst+nbitsd-1
236 endif
237 endif
238 !print *,'SDp ',ival1,ival2,minsd,nbitsd
239 endif ! end of spatial diff section
240 !
241 ! Expand non-missing data values to original grid.
242 !
243 miss1=minval(jfld(1:nonmiss))-1
244 miss2=miss1-1
245 n=0
246 do j=1,ndpts
247 if ( ifldmiss(j).eq.0 ) then
248 n=n+1
249 ifld(j)=jfld(n)
250 elseif ( ifldmiss(j).eq.1 ) then
251 ifld(j)=miss1
252 elseif ( ifldmiss(j).eq.2 ) then
253 ifld(j)=miss2
254 endif
255 enddo
256 !
257 ! Determine Groups to be used.
258 !
259 if ( simple_alg ) then
260 ! set group length to 10 : calculate number of groups
261 ! and length of last group
262 ngroups=ndpts/10
263 glen(1:ngroups)=10
264 itemp=mod(ndpts,10)
265 if (itemp.ne.0) then
266 ngroups=ngroups+1
267 glen(ngroups)=itemp
268 endif
269 else
270 ! Use Dr. Glahn's algorithm for determining grouping.
271 !
272 kfildo=6
273 minpk=10
274 inc=1
275 maxgrps=(ndpts/minpk)+1
276 allocate(jmin(maxgrps))
277 allocate(jmax(maxgrps))
278 allocate(lbit(maxgrps))
279 call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2,
280 & jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit,
281 & kbit,novref,lbitref,ier)
282 !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
283 do ng=1,ngroups
284 glen(ng)=glen(ng)+novref
285 enddo
286 deallocate(jmin)
287 deallocate(jmax)
288 deallocate(lbit)
289 endif
290 !
291 ! For each group, find the group's reference value (min)
292 ! and the number of bits needed to hold the remaining values
293 !
294 n=1
295 do ng=1,ngroups
296 ! how many of each type?
297 num0=count(ifldmiss(n:n+glen(ng)-1) .EQ. 0)
298 num1=count(ifldmiss(n:n+glen(ng)-1) .EQ. 1)
299 num2=count(ifldmiss(n:n+glen(ng)-1) .EQ. 2)
300 if ( num0.eq.0 ) then ! all missing values
301 if ( num1.eq.0 ) then ! all secondary missing
302 gref(ng)=-2
303 gwidth(ng)=0
304 elseif ( num2.eq.0 ) then ! all primary missing
305 gref(ng)=-1
306 gwidth(ng)=0
307 else ! both primary and secondary
308 gref(ng)=0
309 gwidth(ng)=1
310 endif
311 else ! contains some non-missing data
312 ! find max and min values of group
313 gref(ng)=huge(n)
314 imax=-1*huge(n)
315 j=n
316 do lg=1,glen(ng)
317 if ( ifldmiss(j).eq.0 ) then
318 if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j)
319 if (ifld(j).gt.imax) imax=ifld(j)
320 endif
321 j=j+1
322 enddo
323 if (missopt.eq.1) imax=imax+1
324 if (missopt.eq.2) imax=imax+2
325 ! calc num of bits needed to hold data
326 if ( gref(ng).ne.imax ) then
327 temp=alog(real(imax-gref(ng)+1))/alog2
328 gwidth(ng)=ceiling(temp)
329 else
330 gwidth(ng)=0
331 endif
332 endif
333 ! Subtract min from data
334 j=n
335 mtemp=2**gwidth(ng)
336 do lg=1,glen(ng)
337 if (ifldmiss(j).eq.0) then ! non-missing
338 ifld(j)=ifld(j)-gref(ng)
339 elseif (ifldmiss(j).eq.1) then ! primary missing
340 ifld(j)=mtemp-1
341 elseif (ifldmiss(j).eq.2) then ! secondary missing
342 ifld(j)=mtemp-2
343 endif
344 j=j+1
345 enddo
346 ! increment fld array counter
347 n=n+glen(ng)
348 enddo
349 !
350 ! Find max of the group references and calc num of bits needed
351 ! to pack each groups reference value, then
352 ! pack up group reference values
353 !
354 !write(77,*)'GREFS: ',(gref(j),j=1,ngroups)
355 igmax=maxval(gref(1:ngroups))
356 if (missopt.eq.1) igmax=igmax+1
357 if (missopt.eq.2) igmax=igmax+2
358 if (igmax.ne.0) then
359 temp=alog(real(igmax+1))/alog2
360 nbitsgref=ceiling(temp)
361 ! restet the ref values of any "missing only" groups.
362 mtemp=2**nbitsgref
363 do j=1,ngroups
364 if (gref(j).eq.-1) gref(j)=mtemp-1
365 if (gref(j).eq.-2) gref(j)=mtemp-2
366 enddo
367 call g2lib_sbytes(cpack,gref,iofst,nbitsgref,0,ngroups)
368 itemp=nbitsgref*ngroups
369 iofst=iofst+itemp
370 ! Pad last octet with Zeros, if necessary,
371 if (mod(itemp,8).ne.0) then
372 left=8-mod(itemp,8)
373 call g2lib_sbyte(cpack,zero,iofst,left)
374 iofst=iofst+left
375 endif
376 else
377 nbitsgref=0
378 endif
379 !
380 ! Find max/min of the group widths and calc num of bits needed
381 ! to pack each groups width value, then
382 ! pack up group width values
383 !
384 !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
385 iwmax=maxval(gwidth(1:ngroups))
386 ngwidthref=minval(gwidth(1:ngroups))
387 if (iwmax.ne.ngwidthref) then
388 temp=alog(real(iwmax-ngwidthref+1))/alog2
389 nbitsgwidth=ceiling(temp)
390 do i=1,ngroups
391 gwidth(i)=gwidth(i)-ngwidthref
392 enddo
393 call g2lib_sbytes(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
394 itemp=nbitsgwidth*ngroups
395 iofst=iofst+itemp
396 ! Pad last octet with Zeros, if necessary,
397 if (mod(itemp,8).ne.0) then
398 left=8-mod(itemp,8)
399 call g2lib_sbyte(cpack,zero,iofst,left)
400 iofst=iofst+left
401 endif
402 else
403 nbitsgwidth=0
404 gwidth(1:ngroups)=0
405 endif
406 !
407 ! Find max/min of the group lengths and calc num of bits needed
408 ! to pack each groups length value, then
409 ! pack up group length values
410 !
411 !write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
412 ilmax=maxval(glen(1:ngroups-1))
413 nglenref=minval(glen(1:ngroups-1))
414 nglenlast=glen(ngroups)
415 if (ilmax.ne.nglenref) then
416 temp=alog(real(ilmax-nglenref+1))/alog2
417 nbitsglen=ceiling(temp)
418 do i=1,ngroups-1
419 glen(i)=glen(i)-nglenref
420 enddo
421 call g2lib_sbytes(cpack,glen,iofst,nbitsglen,0,ngroups)
422 itemp=nbitsglen*ngroups
423 iofst=iofst+itemp
424 ! Pad last octet with Zeros, if necessary,
425 if (mod(itemp,8).ne.0) then
426 left=8-mod(itemp,8)
427 call g2lib_sbyte(cpack,zero,iofst,left)
428 iofst=iofst+left
429 endif
430 else
431 nbitsglen=0
432 glen(1:ngroups)=0
433 endif
434 !
435 ! For each group, pack data values
436 !
437 !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
438 n=1
439 ij=0
440 do ng=1,ngroups
441 glength=glen(ng)+nglenref
442 if (ng.eq.ngroups ) glength=nglenlast
443 grpwidth=gwidth(ng)+ngwidthref
444 !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
445 if ( grpwidth.ne.0 ) then
446 call g2lib_sbytes(cpack,ifld(n),iofst,grpwidth,0,glength)
447 iofst=iofst+(grpwidth*glength)
448 endif
449 do kk=1,glength
450 ij=ij+1
451 !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
452 enddo
453 n=n+glength
454 enddo
455 ! Pad last octet with Zeros, if necessary,
456 if (mod(iofst,8).ne.0) then
457 left=8-mod(iofst,8)
458 call g2lib_sbyte(cpack,zero,iofst,left)
459 iofst=iofst+left
460 endif
461 lcpack=iofst/8
462 !
463 if ( allocated(ifld) ) deallocate(ifld)
464 if ( allocated(jfld) ) deallocate(jfld)
465 if ( allocated(ifldmiss) ) deallocate(ifldmiss)
466 if ( allocated(gref) ) deallocate(gref)
467 if ( allocated(gwidth) ) deallocate(gwidth)
468 if ( allocated(glen) ) deallocate(glen)
469 !else ! Constant field ( max = min )
470 ! nbits=0
471 ! lcpack=0
472 ! nbitsgref=0
473 ! ngroups=0
474 !endif
476 !
477 ! Fill in ref value and number of bits in Template 5.2
478 !
479 call mkieee(rmin,ref,1) ! ensure reference value is IEEE format
480 ! call g2lib_gbyte(ref,idrstmpl(1),0,32)
481 iref=transfer(ref,iref)
482 idrstmpl(1)=iref
483 idrstmpl(4)=nbitsgref
484 idrstmpl(5)=0 ! original data were reals
485 idrstmpl(6)=1 ! general group splitting
486 idrstmpl(10)=ngroups ! Number of groups
487 idrstmpl(11)=ngwidthref ! reference for group widths
488 idrstmpl(12)=nbitsgwidth ! num bits used for group widths
489 idrstmpl(13)=nglenref ! Reference for group lengths
490 idrstmpl(14)=1 ! length increment for group lengths
491 idrstmpl(15)=nglenlast ! True length of last group
492 idrstmpl(16)=nbitsglen ! num bits used for group lengths
493 if (idrsnum.eq.3) then
494 idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial
495 ! differencing values
496 endif
498 return
499 end