ungrib/src/ngl/g2/misspack.f

   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 !$$$
  60
  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
  66
  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.
  75
  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 sbyte(cpack,ival1,iofst,nbitsd)
 208                  iofst=iofst+nbitsd
 209               else
 210                  call sbyte(cpack,1,iofst,1)
 211                  iofst=iofst+1
 212                  call 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 sbyte(cpack,ival2,iofst,nbitsd)
 219                     iofst=iofst+nbitsd
 220                  else
 221                     call sbyte(cpack,1,iofst,1)
 222                     iofst=iofst+1
 223                     call 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 sbyte(cpack,minsd,iofst,nbitsd)
 230                  iofst=iofst+nbitsd
 231               else
 232                  call sbyte(cpack,1,iofst,1)
 233                  iofst=iofst+1
 234                  call 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 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 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 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 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 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 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 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 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
 475
 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 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
 497
 498       return
 499       end