ungrib build
[WPS.git] / ungrib / src / ngl / g2 / misspack.f
blobb338307fb3a91c5f98cabd7d27c7e46c1ec1e099
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
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:
17 ! 2000-06-21 Gilbert
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
26 ! Must equal 2 or 3.
27 ! idrstmpl - Contains the array of values for Data Representation
28 ! Template 5.2 or 5.3
29 ! (1) = Reference value - ignored on input
30 ! (2) = Binary Scale Factor
31 ! (3) = Decimal Scale Factor
32 ! .
33 ! .
34 ! (7) = Missing value management
35 ! (8) = Primary missing value
36 ! (9) = Secondary missing value
37 ! .
38 ! .
39 ! (17) = Order of Spatial Differencing ( 1 or 2 )
40 ! .
41 ! .
43 ! OUTPUT ARGUMENT LIST:
44 ! idrstmpl - Contains the array of values for Data Representation
45 ! Template 5.3
46 ! (1) = Reference value - set by misspack routine.
47 ! (2) = Binary Scale Factor - unchanged from input
48 ! (3) = Decimal Scale Factor - unchanged from input
49 ! .
50 ! .
51 ! cpack - The packed data field (character*1 array)
52 ! lcpack - length of packed field cpack().
54 ! REMARKS: None
56 ! ATTRIBUTES:
57 ! LANGUAGE: XL Fortran 90
58 ! MACHINE: IBM SP
60 !$$$
61 use intmath
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
68 real(4) :: ref, rmin4
69 integer(4) :: iref
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
75 logical :: simple_alg
76 logical :: have_rmin
78 simple_alg = .false.
79 have_rmin = .false.
80 bscale=2.0**real(-idrstmpl(2))
81 dscale=10.0**real(idrstmpl(3))
82 missopt=idrstmpl(7)
83 if ( missopt.ne.1 .AND. missopt.ne.2 ) then
84 print *,'misspack: Unrecognized option.'
85 lcpack=-1
86 return
87 else ! Get missing values
88 call rdieee(idrstmpl(8),rmissp,1)
89 if (missopt.eq.2) call rdieee(idrstmpl(9),rmisss,1)
90 endif
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))
96 c rmin=huge(rmin)
98 if ( missopt .eq. 1 ) then ! Primary missing value only
99 do j=1,ndpts
100 if (fld(j).eq.rmissp) then
101 ifldmiss(j)=1
102 else
103 ifldmiss(j)=0
104 if(have_rmin) then
105 if (fld(j).lt.rmin) rmin=fld(j)
106 else
107 rmin=fld(j)
108 have_rmin=.true.
109 endif
110 endif
111 enddo
112 if(.not.have_rmin) rmin=rmissp
113 endif
114 if ( missopt .eq. 2 ) then ! Primary and secondary missing values
115 do j=1,ndpts
116 if (fld(j).eq.rmissp) then
117 ifldmiss(j)=1
118 elseif (fld(j).eq.rmisss) then
119 ifldmiss(j)=2
120 else
121 ifldmiss(j)=0
122 if(have_rmin) then
123 if (fld(j).lt.rmin) rmin=fld(j)
124 else
125 rmin=fld(j)
126 have_rmin=.true.
127 endif
128 endif
129 if(.not.have_rmin) rmin=rmissp
130 enddo
131 endif
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
141 iofst=0
142 allocate(ifld(ndpts))
143 allocate(jfld(ndpts))
144 allocate(gref(ndpts))
145 allocate(gwidth(ndpts))
146 allocate(glen(ndpts))
148 ! Scale original data
150 nonmiss=0
151 if (idrstmpl(2).eq.0) then ! No binary scaling
152 imin=nint(rmin*dscale)
153 !imax=nint(rmax*dscale)
154 rmin=real(imin)
155 do j=1,ndpts
156 if (ifldmiss(j).eq.0) then
157 nonmiss=nonmiss+1
158 jfld(nonmiss)=max(0,nint(fld(j)*dscale)-imin)
159 endif
160 enddo
161 else ! Use binary scaling factor
162 rmin=rmin*dscale
163 !rmax=rmax*dscale
164 do j=1,ndpts
165 if (ifldmiss(j).eq.0) then
166 nonmiss=nonmiss+1
167 jfld(nonmiss)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
168 endif
169 enddo
170 endif
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
177 if(nonmiss<1) then
178 ival1=1.0
179 else
180 ival1=jfld(1)
181 endif
182 do j=nonmiss,2,-1
183 jfld(j)=jfld(j)-jfld(j-1)
184 enddo
185 if(nonmiss>0) jfld(1)=0
186 elseif (idrstmpl(17).eq.2) then ! second order
187 if(nonmiss==1) then
188 ival1=jfld(1)
189 ival2=jfld(1)
190 elseif(nonmiss<1) then
191 ival1=1.0
192 ival2=1.0
193 else
194 ival1=jfld(1)
195 ival2=jfld(2)
196 endif
197 do j=nonmiss,3,-1
198 jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2)
199 enddo
200 if(nonmiss>=1) jfld(1)=0
201 if(nonmiss>=2) jfld(2)=0
202 endif
204 ! subtract min value from spatial diff field
206 isd=idrstmpl(17)+1
207 minsd=minval(jfld(isd:nonmiss))
208 do j=isd,nonmiss
209 jfld(j)=jfld(j)-minsd
210 enddo
212 ! find num of bits need to store minsd and add 1 extra bit
213 ! to indicate sign
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 )
221 maxorig=ival1
222 if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2
223 temp=i1log2(maxorig)
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
230 ! data section.
232 if (nbitsd.ne.0) then
233 ! pack first original value
234 if (ival1.ge.0) then
235 call sbyte(cpack,ival1,iofst,nbitsd)
236 iofst=iofst+nbitsd
237 else
238 call sbyte(cpack,1,iofst,1)
239 iofst=iofst+1
240 call sbyte(cpack,iabs(ival1),iofst,nbitsd-1)
241 iofst=iofst+nbitsd-1
242 endif
243 if (idrstmpl(17).eq.2) then
244 ! pack second original value
245 if (ival2.ge.0) then
246 call sbyte(cpack,ival2,iofst,nbitsd)
247 iofst=iofst+nbitsd
248 else
249 call sbyte(cpack,1,iofst,1)
250 iofst=iofst+1
251 call sbyte(cpack,iabs(ival2),iofst,nbitsd-1)
252 iofst=iofst+nbitsd-1
253 endif
254 endif
255 ! pack overall min of spatial differences
256 if (minsd.ge.0) then
257 call sbyte(cpack,minsd,iofst,nbitsd)
258 iofst=iofst+nbitsd
259 else
260 call sbyte(cpack,1,iofst,1)
261 iofst=iofst+1
262 call sbyte(cpack,iabs(minsd),iofst,nbitsd-1)
263 iofst=iofst+nbitsd-1
264 endif
265 endif
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
272 miss2=miss1-1
274 do j=1,ndpts
275 if ( ifldmiss(j).eq.0 ) then
276 n=n+1
277 ifld(j)=jfld(n)
278 elseif ( ifldmiss(j).eq.1 ) then
279 ifld(j)=miss1
280 elseif ( ifldmiss(j).eq.2 ) then
281 ifld(j)=miss2
282 endif
283 enddo
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
291 ngroups=ndpts/10
292 glen(1:ngroups)=10
293 itemp=mod(ndpts,10)
294 if (itemp.ne.0) then
295 ngroups=ngroups+1
296 glen(ngroups)=itemp
297 endif
298 else
299 ! Use Dr. Glahn's algorithm for determining grouping.
301 kfildo=6
302 minpk=10
303 inc=1
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
312 do ng=1,ngroups
313 glen(ng)=glen(ng)+novref
314 enddo
315 deallocate(jmin)
316 deallocate(jmax)
317 deallocate(lbit)
318 endif
320 ! For each group, find the group's reference value (min)
321 ! and the number of bits needed to hold the remaining values
324 do ng=1,ngroups
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
331 gref(ng)=-2
332 gwidth(ng)=0
333 elseif ( num2.eq.0 ) then ! all primary missing
334 gref(ng)=-1
335 gwidth(ng)=0
336 else ! both primary and secondary
337 gref(ng)=0
338 gwidth(ng)=1
339 endif
340 else ! contains some non-missing data
341 ! find max and min values of group
342 gref(ng)=huge(n)
343 imax=-1*huge(n)
345 do lg=1,glen(ng)
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)
349 endif
350 j=j+1
351 enddo
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)
358 else
359 gwidth(ng)=0
360 endif
361 endif
362 ! Subtract min from data
364 mtemp=2**gwidth(ng)
365 do lg=1,glen(ng)
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
369 ifld(j)=mtemp-1
370 elseif (ifldmiss(j).eq.2) then ! secondary missing
371 ifld(j)=mtemp-2
372 endif
373 j=j+1
374 enddo
375 ! increment fld array counter
376 n=n+glen(ng)
377 enddo
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
387 if (igmax.ne.0) then
388 temp=i1log2(igmax)
389 nbitsgref=ceiling(temp)
390 ! restet the ref values of any "missing only" groups.
391 mtemp=2**nbitsgref
392 do j=1,ngroups
393 if (gref(j).eq.-1) gref(j)=mtemp-1
394 if (gref(j).eq.-2) gref(j)=mtemp-2
395 enddo
396 call sbytes(cpack,gref,iofst,nbitsgref,0,ngroups)
397 itemp=nbitsgref*ngroups
398 iofst=iofst+itemp
399 ! Pad last octet with Zeros, if necessary,
400 if (mod(itemp,8).ne.0) then
401 left=8-mod(itemp,8)
402 call sbyte(cpack,zero,iofst,left)
403 iofst=iofst+left
404 endif
405 else
406 nbitsgref=0
407 endif
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)
419 do i=1,ngroups
420 gwidth(i)=gwidth(i)-ngwidthref
421 enddo
422 call sbytes(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
423 itemp=nbitsgwidth*ngroups
424 iofst=iofst+itemp
425 ! Pad last octet with Zeros, if necessary,
426 if (mod(itemp,8).ne.0) then
427 left=8-mod(itemp,8)
428 call sbyte(cpack,zero,iofst,left)
429 iofst=iofst+left
430 endif
431 else
432 nbitsgwidth=0
433 gwidth(1:ngroups)=0
434 endif
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))
443 if(ngroups>0) then
444 nglenlast=glen(ngroups)
445 else
446 nglenlast=0
447 endif
448 if (ilmax.ne.nglenref) then
449 temp=i1log2(ilmax-nglenref)
450 nbitsglen=ceiling(temp)
451 do i=1,ngroups-1
452 glen(i)=glen(i)-nglenref
453 enddo
454 call sbytes(cpack,glen,iofst,nbitsglen,0,ngroups)
455 itemp=nbitsglen*ngroups
456 iofst=iofst+itemp
457 ! Pad last octet with Zeros, if necessary,
458 if (mod(itemp,8).ne.0) then
459 left=8-mod(itemp,8)
460 call sbyte(cpack,zero,iofst,left)
461 iofst=iofst+left
462 endif
463 else
464 nbitsglen=0
465 glen(1:ngroups)=0
466 endif
468 ! For each group, pack data values
470 !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
472 ij=0
473 do ng=1,ngroups
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)
481 endif
482 do kk=1,glength
483 ij=ij+1
484 !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
485 enddo
486 n=n+glength
487 enddo
488 ! Pad last octet with Zeros, if necessary,
489 if (mod(iofst,8).ne.0) then
490 left=8-mod(iofst,8)
491 call sbyte(cpack,zero,iofst,left)
492 iofst=iofst+left
493 endif
494 lcpack=iofst/8
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 )
503 ! nbits=0
504 ! lcpack=0
505 ! nbitsgref=0
506 ! ngroups=0
507 !endif
510 ! Fill in ref value and number of bits in Template 5.2
512 rmin4 = rmin
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)
516 idrstmpl(1)=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
530 endif
532 return