ungrib/src/ngl/g2/compack.f

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