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