1 subroutine misspack
(fld
,ndpts
,idrsnum
,idrstmpl
,cpack
,lcpack
)
2 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
5 ! PRGMMR
: Gilbert ORG
: W
/NP11
DATE: 2000-06-21
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
.
16 ! PROGRAM HISTORY LOG
:
18 ! 2004-12-29 Gilbert
- Corrected bug when encoding secondary missing values
.
19 ! 2012-05-10 Boi Vuong Added variable rmin4
for 4 byte
real
21 ! USAGE
: CALL misspack
(fld
,ndpts
,idrsnum
,idrstmpl
,cpack
,lcpack
)
22 ! INPUT ARGUMENT LIST
:
23 ! fld
() - Contains the data values
to pack
24 ! ndpts
- The number of data values in array fld
()
25 ! idrsnum
- Data Representation Template number
5.N
27 ! idrstmpl
- Contains the array of values
for Data Representation
29 ! (1) = Reference value
- ignored on input
30 ! (2) = Binary Scale Factor
31 ! (3) = Decimal Scale Factor
34 ! (7) = Missing value management
35 ! (8) = Primary missing value
36 ! (9) = Secondary missing value
39 ! (17) = Order of Spatial Differencing
( 1 or
2 )
43 ! OUTPUT ARGUMENT LIST
:
44 ! idrstmpl
- Contains the array of values
for Data Representation
46 ! (1) = Reference value
- set by misspack routine
.
47 ! (2) = Binary Scale Factor
- unchanged from input
48 ! (3) = Decimal Scale Factor
- unchanged from input
51 ! cpack
- The packed data field
(character*1 array
)
52 ! lcpack
- length of packed field cpack
().
57 ! LANGUAGE
: XL Fortran
90
62 integer,intent
(in
) :: ndpts
,idrsnum
63 real,intent
(in
) :: fld
(ndpts
)
64 character(len
=1),intent
(out
) :: cpack
(*)
65 integer,intent
(inout
) :: idrstmpl
(*)
66 integer,intent
(out
) :: lcpack
70 integer,allocatable
:: ifld
(:),ifldmiss
(:),jfld
(:)
71 integer,allocatable
:: jmin
(:),jmax
(:),lbit
(:)
72 integer,parameter :: zero
=0
73 integer,allocatable
:: gref
(:),gwidth
(:),glen
(:)
74 integer :: glength
,grpwidth
80 bscale
=2.0**real(-idrstmpl
(2))
81 dscale
=10.0**real(idrstmpl
(3))
83 if ( missopt
.ne
.1 .AND
. missopt
.ne
.2 ) then
84 print
*,'misspack: Unrecognized option.'
87 else ! Get missing values
88 call rdieee
(idrstmpl
(8),rmissp
,1)
89 if (missopt
.eq
.2) call rdieee
(idrstmpl
(9),rmisss
,1)
92 ! Find min value of non
-missing values in the data
,
93 ! AND set up missing value mapping of the field
.
95 allocate
(ifldmiss
(ndpts
))
98 if ( missopt
.eq
. 1 ) then ! Primary missing value only
100 if (fld
(j
).eq
.rmissp
) then
105 if (fld
(j
).lt
.rmin
) rmin
=fld
(j
)
112 if(.not
.have_rmin
) rmin
=rmissp
114 if ( missopt
.eq
. 2 ) then ! Primary and secondary missing values
116 if (fld
(j
).eq
.rmissp
) then
118 elseif
(fld
(j
).eq
.rmisss
) then
123 if (fld
(j
).lt
.rmin
) rmin
=fld
(j
)
129 if(.not
.have_rmin
) rmin
=rmissp
133 ! Allocate work arrays
:
134 ! Note
: -ifldmiss
(j
),j
=1,ndpts is a map of original field indicating
135 ! which of the original data values
136 ! are primary missing
(1), sencondary missing
(2) or non
-missing
(0).
137 ! -jfld
(j
),j
=1,nonmiss is a subarray of just the non
-missing values from
138 ! the original field
.
140 !if (rmin
.ne
.rmax
) then
142 allocate
(ifld
(ndpts
))
143 allocate
(jfld
(ndpts
))
144 allocate
(gref
(ndpts
))
145 allocate
(gwidth
(ndpts
))
146 allocate
(glen
(ndpts
))
148 ! Scale original data
151 if (idrstmpl
(2).eq
.0) then ! No binary scaling
152 imin
=nint
(rmin*dscale
)
153 !imax
=nint
(rmax*dscale
)
156 if (ifldmiss
(j
).eq
.0) then
158 jfld
(nonmiss
)=max
(0,nint
(fld
(j
)*dscale
)-imin
)
161 else ! Use binary scaling factor
165 if (ifldmiss
(j
).eq
.0) then
167 jfld
(nonmiss
)=max
(0,nint
(((fld
(j
)*dscale
)-rmin
)*bscale
))
172 ! Calculate Spatial differences
, if using DRS Template
5.3
174 if (idrsnum
.eq
.3) then ! spatial differences
175 if (idrstmpl
(17).ne
.1.and
.idrstmpl
(17).ne
.2) idrstmpl
(17)=2
176 if (idrstmpl
(17).eq
.1) then ! first order
183 jfld
(j
)=jfld
(j
)-jfld
(j
-1)
185 if(nonmiss
>0) jfld
(1)=0
186 elseif
(idrstmpl
(17).eq
.2) then ! second order
190 elseif
(nonmiss
<1) then
198 jfld
(j
)=jfld
(j
)-(2*jfld
(j
-1))+jfld
(j
-2)
200 if(nonmiss
>=1) jfld
(1)=0
201 if(nonmiss
>=2) jfld
(2)=0
204 ! subtract min value from spatial diff field
207 minsd
=minval
(jfld
(isd
:nonmiss
))
209 jfld
(j
)=jfld
(j
)-minsd
212 ! find num of bits need
to store minsd and add
1 extra bit
215 temp
=i1log2
(abs
(minsd
))
216 nbitsd
=ceiling
(temp
)+1
218 ! find num of bits need
to store ifld
(1) ( and ifld
(2)
219 ! if using
2nd order differencing
)
222 if (idrstmpl
(17).eq
.2.and
.ival2
.gt
.ival1
) maxorig
=ival2
224 nbitorig
=ceiling
(temp
)+1
225 if (nbitorig
.gt
.nbitsd
) nbitsd
=nbitorig
226 ! increase number of bits
to even multiple of
8 ( octet
)
227 if (mod
(nbitsd
,8).ne
.0) nbitsd
=nbitsd
+(8-mod
(nbitsd
,8))
229 ! Store extra spatial differencing info into the packed
232 if (nbitsd
.ne
.0) then
233 ! pack first original value
235 call sbyte
(cpack
,ival1
,iofst
,nbitsd
)
238 call sbyte
(cpack
,1,iofst
,1)
240 call sbyte
(cpack
,iabs
(ival1
),iofst
,nbitsd
-1)
243 if (idrstmpl
(17).eq
.2) then
244 ! pack second original value
246 call sbyte
(cpack
,ival2
,iofst
,nbitsd
)
249 call sbyte
(cpack
,1,iofst
,1)
251 call sbyte
(cpack
,iabs
(ival2
),iofst
,nbitsd
-1)
255 ! pack overall min of spatial differences
257 call sbyte
(cpack
,minsd
,iofst
,nbitsd
)
260 call sbyte
(cpack
,1,iofst
,1)
262 call sbyte
(cpack
,iabs
(minsd
),iofst
,nbitsd
-1)
266 !print
*,'SDp ',ival1
,ival2
,minsd
,nbitsd
267 endif ! end of spatial diff section
269 ! Expand non
-missing data values
to original grid
.
271 miss1
=minval
(jfld
(1:nonmiss
))-1
275 if ( ifldmiss
(j
).eq
.0 ) then
278 elseif
( ifldmiss
(j
).eq
.1 ) then
280 elseif
( ifldmiss
(j
).eq
.2 ) then
284 if(ndpts
<2) simple_alg
=.true
.
286 ! Determine Groups
to be used
.
288 if ( simple_alg
) then
289 ! set group length
to 10 : calculate number of groups
290 ! and length of last group
299 ! Use Dr
. Glahn
's algorithm for determining grouping.
304 maxgrps=(ndpts/minpk)+1
305 allocate(jmin(maxgrps))
306 allocate(jmax(maxgrps))
307 allocate(lbit(maxgrps))
308 call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2,
309 & jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit,
310 & kbit,novref,lbitref,ier)
311 !print *,'SAGier
= ',ier,ibit,jbit,kbit,novref,lbitref
313 glen(ng)=glen(ng)+novref
320 ! For each group, find the group's reference value
(min
)
321 ! and the number of bits needed
to hold the remaining values
325 ! how many of each type?
326 num0
=count
(ifldmiss
(n
:n
+glen
(ng
)-1) .EQ
. 0)
327 num1
=count
(ifldmiss
(n
:n
+glen
(ng
)-1) .EQ
. 1)
328 num2
=count
(ifldmiss
(n
:n
+glen
(ng
)-1) .EQ
. 2)
329 if ( num0
.eq
.0 ) then ! all missing values
330 if ( num1
.eq
.0 ) then ! all secondary missing
333 elseif
( num2
.eq
.0 ) then ! all primary missing
336 else ! both primary and secondary
340 else ! contains some non
-missing data
341 ! find max and min values of group
346 if ( ifldmiss
(j
).eq
.0 ) then
347 if (ifld
(j
).lt
.gref
(ng
)) gref
(ng
)=ifld
(j
)
348 if (ifld
(j
).gt
.imax
) imax
=ifld
(j
)
352 if (missopt
.eq
.1) imax
=imax
+1
353 if (missopt
.eq
.2) imax
=imax
+2
354 ! calc num of bits needed
to hold data
355 if ( gref
(ng
).ne
.imax
) then
356 temp
=i1log2
(imax
-gref
(ng
))
357 gwidth
(ng
)=ceiling
(temp
)
362 ! Subtract min from data
366 if (ifldmiss
(j
).eq
.0) then ! non
-missing
367 ifld
(j
)=ifld
(j
)-gref
(ng
)
368 elseif
(ifldmiss
(j
).eq
.1) then ! primary missing
370 elseif
(ifldmiss
(j
).eq
.2) then ! secondary missing
375 ! increment fld array counter
379 ! Find max of the group references and calc num of bits needed
380 ! to pack each groups reference value
, then
381 ! pack up group reference values
383 !write(77,*)'GREFS: ',(gref
(j
),j
=1,ngroups
)
384 igmax
=maxval
(gref
(1:ngroups
))
385 if (missopt
.eq
.1) igmax
=igmax
+1
386 if (missopt
.eq
.2) igmax
=igmax
+2
389 nbitsgref
=ceiling
(temp
)
390 ! restet the ref values of any
"missing only" groups
.
393 if (gref
(j
).eq
.-1) gref
(j
)=mtemp
-1
394 if (gref
(j
).eq
.-2) gref
(j
)=mtemp
-2
396 call sbytes
(cpack
,gref
,iofst
,nbitsgref
,0,ngroups
)
397 itemp
=nbitsgref*ngroups
399 ! Pad last octet with Zeros
, if necessary
,
400 if (mod
(itemp
,8).ne
.0) then
402 call sbyte
(cpack
,zero
,iofst
,left
)
409 ! Find max
/min of the group widths and calc num of bits needed
410 ! to pack each groups width value
, then
411 ! pack up group width values
413 !write(77,*)'GWIDTHS: ',(gwidth
(j
),j
=1,ngroups
)
414 iwmax
=maxval
(gwidth
(1:ngroups
))
415 ngwidthref
=minval
(gwidth
(1:ngroups
))
416 if (iwmax
.ne
.ngwidthref
) then
417 temp
=i1log2
(iwmax
-ngwidthref
)
418 nbitsgwidth
=ceiling
(temp
)
420 gwidth
(i
)=gwidth
(i
)-ngwidthref
422 call sbytes
(cpack
,gwidth
,iofst
,nbitsgwidth
,0,ngroups
)
423 itemp
=nbitsgwidth*ngroups
425 ! Pad last octet with Zeros
, if necessary
,
426 if (mod
(itemp
,8).ne
.0) then
428 call sbyte
(cpack
,zero
,iofst
,left
)
436 ! Find max
/min of the group lengths and calc num of bits needed
437 ! to pack each groups length value
, then
438 ! pack up group length values
440 !write(77,*)'GLENS: ',(glen
(j
),j
=1,ngroups
)
441 ilmax
=maxval
(glen
(1:ngroups
-1))
442 nglenref
=minval
(glen
(1:ngroups
-1))
444 nglenlast
=glen
(ngroups
)
448 if (ilmax
.ne
.nglenref
) then
449 temp
=i1log2
(ilmax
-nglenref
)
450 nbitsglen
=ceiling
(temp
)
452 glen
(i
)=glen
(i
)-nglenref
454 call sbytes
(cpack
,glen
,iofst
,nbitsglen
,0,ngroups
)
455 itemp
=nbitsglen*ngroups
457 ! Pad last octet with Zeros
, if necessary
,
458 if (mod
(itemp
,8).ne
.0) then
460 call sbyte
(cpack
,zero
,iofst
,left
)
468 ! For each group
, pack data values
470 !write(77,*)'IFLDS: ',(ifld
(j
),j
=1,ndpts
)
474 glength
=glen
(ng
)+nglenref
475 if (ng
.eq
.ngroups
) glength
=nglenlast
476 grpwidth
=gwidth
(ng
)+ngwidthref
477 !write(77,*)'NGP ',ng
,grpwidth
,glength
,gref
(ng
)
478 if ( grpwidth
.ne
.0 ) then
479 call sbytes
(cpack
,ifld
(n
),iofst
,grpwidth
,0,glength
)
480 iofst
=iofst
+(grpwidth*glength
)
484 !write(77,*)'SAG ',ij
,fld
(ij
),ifld
(ij
),gref
(ng
),bscale
,rmin
,dscale
488 ! Pad last octet with Zeros
, if necessary
,
489 if (mod
(iofst
,8).ne
.0) then
491 call sbyte
(cpack
,zero
,iofst
,left
)
496 if ( allocated
(ifld
) ) deallocate
(ifld
)
497 if ( allocated
(jfld
) ) deallocate
(jfld
)
498 if ( allocated
(ifldmiss
) ) deallocate
(ifldmiss
)
499 if ( allocated
(gref
) ) deallocate
(gref
)
500 if ( allocated
(gwidth
) ) deallocate
(gwidth
)
501 if ( allocated
(glen
) ) deallocate
(glen
)
502 !else ! Constant field
( max
= min
)
510 ! Fill in ref value and number of bits in Template
5.2
513 call mkieee
(rmin4
,ref
,1) ! ensure reference value is IEEE
format
514 ! call gbyte
(ref
,idrstmpl
(1),0,32)
515 iref
=transfer
(ref
,iref
)
517 idrstmpl
(4)=nbitsgref
518 idrstmpl
(5)=0 ! original data were reals
519 idrstmpl
(6)=1 ! general group splitting
520 idrstmpl
(10)=ngroups
! Number of groups
521 idrstmpl
(11)=ngwidthref
! reference
for group widths
522 idrstmpl
(12)=nbitsgwidth
! num bits used
for group widths
523 idrstmpl
(13)=nglenref
! Reference
for group lengths
524 idrstmpl
(14)=1 ! length increment
for group lengths
525 idrstmpl
(15)=nglenlast
! True length of last group
526 idrstmpl
(16)=nbitsglen
! num bits used
for group lengths
527 if (idrsnum
.eq
.3) then
528 idrstmpl
(18)=nbitsd
/8 ! num bits used
for extra spatial
529 ! differencing values