ungrib/src/ngl/g2/comunpack.f

   1       subroutine comunpack(cpack,len,lensec,idrsnum,idrstmpl,ndpts,
   2      &                     fld,ier)
   3 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   4 !                .      .    .                                       .
   5 ! SUBPROGRAM:    comunpack
   6 !   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2000-06-21
   7 !
   8 ! ABSTRACT: This subroutine unpacks a data field that was packed using a
   9 !   complex packing algorithm as defined in the GRIB2 documention,
  10 !   using info from the GRIB2 Data Representation Template 5.2 or 5.3.
  11 !   Supports GRIB2 complex packing templates with or without
  12 !   spatial differences (i.e. DRTs 5.2 and 5.3).
  13 !
  14 ! PROGRAM HISTORY LOG:
  15 ! 2000-06-21  Gilbert
  16 ! 2004-12-29  Gilbert  -  Added test ( provided by Arthur Taylor/MDL )
  17 !                         to verify that group widths and lengths are
  18 !                         consistent with section length.
  19 ! 2016-02-26              update unpacking for template 5.3
  20 !
  21 ! USAGE:    CALL comunpack(cpack,len,lensec,idrsnum,idrstmpl,ndpts,fld,ier)
  22 !   INPUT ARGUMENT LIST:
  23 !     cpack    - The packed data field (character*1 array)
  24 !     len      - length of packed field cpack().
  25 !     lensec   - length of section 7 (used for error checking).
  26 !     idrsnum  - Data Representation Template number 5.N
  27 !                Must equal 2 or 3.
  28 !     idrstmpl - Contains the array of values for Data Representation
  29 !                Template 5.2 or 5.3
  30 !     ndpts    - The number of data values to unpack
  31 !
  32 !   OUTPUT ARGUMENT LIST:
  33 !     fld()    - Contains the unpacked data values
  34 !     ier      - Error return:
  35 !                  0 = OK
  36 !                  1 = Problem - inconsistent group lengths of widths.
  37 !
  38 ! REMARKS: None
  39 !
  40 ! ATTRIBUTES:
  41 !   LANGUAGE: XL Fortran 90
  42 !   MACHINE:  IBM SP
  43 !
  44 !$$$
  45
  46       character(len=1),intent(in) :: cpack(len)
  47       integer,intent(in) :: ndpts,len
  48       integer,intent(in) :: idrstmpl(*)
  49       real,intent(out) :: fld(ndpts)
  50
  51       integer,allocatable :: ifld(:),ifldmiss(:)
  52       integer(4) :: ieee
  53       integer,allocatable :: gref(:),gwidth(:),glen(:)
  54       real :: ref,bscale,dscale,rmiss1,rmiss2
  55 !      real :: fldo(6045)
  56       integer :: totBit, totLen
  57
  58       ier=0
  59       !print *,'IDRSTMPL: ',(idrstmpl(j),j=1,16)
  60       ieee = idrstmpl(1)
  61       call rdieee(ieee,ref,1)
  62       bscale = 2.0**real(idrstmpl(2))
  63       dscale = 10.0**real(-idrstmpl(3))
  64       nbitsgref = idrstmpl(4)
  65       itype = idrstmpl(5)
  66       ngroups = idrstmpl(10)
  67       nbitsgwidth = idrstmpl(12)
  68       nbitsglen = idrstmpl(16)
  69       if (idrsnum.eq.3) then
  70          nbitsd=idrstmpl(18)*8
  71       endif
  72
  73       !   Constant field
  74
  75       if (ngroups.eq.0) then
  76          do j=1,ndpts
  77            fld(j)=ref
  78          enddo
  79          return
  80       endif
  81
  82       iofst=0
  83       allocate(ifld(ndpts),stat=is)
  84       !print *,'ALLOC ifld: ',is,ndpts
  85       allocate(gref(ngroups),stat=is)
  86       !print *,'ALLOC gref: ',is
  87       allocate(gwidth(ngroups),stat=is)
  88       !print *,'ALLOC gwidth: ',is
  89 !
  90 !  Get missing values, if supplied
  91 !
  92       if ( idrstmpl(7).eq.1 ) then
  93          if (itype.eq.0) then
  94             call rdieee(idrstmpl(8),rmiss1,1)
  95          else
  96             rmiss1=real(idrstmpl(8))
  97          endif
  98       elseif ( idrstmpl(7).eq.2 ) then
  99          if (itype.eq.0) then
 100             call rdieee(idrstmpl(8),rmiss1,1)
 101             call rdieee(idrstmpl(9),rmiss2,1)
 102          else
 103             rmiss1=real(idrstmpl(8))
 104             rmiss2=real(idrstmpl(9))
 105          endif
 106       endif
 107       !print *,'RMISSs: ',rmiss1,rmiss2,ref
 108 !
 109 !  Extract Spatial differencing values, if using DRS Template 5.3
 110 !
 111       if (idrsnum.eq.3) then
 112          if (nbitsd.ne.0) then
 113               call gbyte(cpack,ival1,iofst,nbitsd)
 114               iofst=iofst+nbitsd
 115               if (idrstmpl(17).eq.2) then
 116                  call gbyte(cpack,ival2,iofst,nbitsd)
 117                  iofst=iofst+nbitsd
 118               endif
 119               call gbyte(cpack,isign,iofst,1)
 120               iofst=iofst+1
 121               call gbyte(cpack,minsd,iofst,nbitsd-1)
 122               iofst=iofst+nbitsd-1
 123               if (isign.eq.1) minsd=-minsd
 124          else
 125               ival1=0
 126               ival2=0
 127               minsd=0
 128          endif
 129        !print *,'SDu ',ival1,ival2,minsd,nbitsd
 130       endif
 131 !
 132 !  Extract Each Group's reference value
 133 !
 134       !print *,'SAG1: ',nbitsgref,ngroups,iofst
 135       if (nbitsgref.ne.0) then
 136          call gbytes(cpack,gref,iofst,nbitsgref,0,ngroups)
 137          itemp=nbitsgref*ngroups
 138          iofst=iofst+(itemp)
 139          if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
 140       else
 141          gref(1:ngroups)=0
 142       endif
 143       !write(78,*)'GREFs: ',(gref(j),j=1,ngroups)
 144 !
 145 !  Extract Each Group's bit width
 146 !
 147       !print *,'SAG2: ',nbitsgwidth,ngroups,iofst,idrstmpl(11)
 148       if (nbitsgwidth.ne.0) then
 149          call gbytes(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
 150          itemp=nbitsgwidth*ngroups
 151          iofst=iofst+(itemp)
 152          if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
 153       else
 154          gwidth(1:ngroups)=0
 155       endif
 156       do j=1,ngroups
 157         gwidth(j)=gwidth(j)+idrstmpl(11)
 158       enddo
 159       !write(78,*)'GWIDTHs: ',(gwidth(j),j=1,ngroups)
 160 !
 161 !  Extract Each Group's length (number of values in each group)
 162 !
 163       allocate(glen(ngroups),stat=is)
 164       !print *,'ALLOC glen: ',is
 165       !print *,'SAG3: ',nbitsglen,ngroups,iofst,idrstmpl(14),idrstmpl(13)
 166       if (nbitsglen.ne.0) then
 167          call gbytes(cpack,glen,iofst,nbitsglen,0,ngroups)
 168          itemp=nbitsglen*ngroups
 169          iofst=iofst+(itemp)
 170          if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8))
 171       else
 172          glen(1:ngroups)=0
 173       endif
 174       do j=1,ngroups
 175         glen(j)=(glen(j)*idrstmpl(14))+idrstmpl(13)
 176       enddo
 177       glen(ngroups)=idrstmpl(15)
 178       !write(78,*)'GLENs: ',(glen(j),j=1,ngroups)
 179       !print *,'GLENsum: ',sum(glen)
 180 !
 181 !  Test to see if the group widths and lengths are consistent with number of
 182 !  values, and length of section 7.
 183 !
 184       totBit = 0
 185       totLen = 0
 186       do j=1,ngroups
 187         totBit = totBit + (gwidth(j)*glen(j));
 188         totLen = totLen + glen(j);
 189       enddo
 190       if (totLen .NE. ndpts) then
 191         ier=1
 192         return
 193       endif
 194       if ( (totBit/8) .GT. lensec) then
 195         ier=1
 196         return
 197       endif
 198 !
 199 !  For each group, unpack data values
 200 !
 201       if ( idrstmpl(7).eq.0 ) then        ! no missing values
 202          n=1
 203          do j=1,ngroups
 204          !write(78,*)'NGP ',j,gwidth(j),glen(j),gref(j)
 205            if (gwidth(j).ne.0) then
 206              call gbytes(cpack,ifld(n),iofst,gwidth(j),0,glen(j))
 207              do k=1,glen(j)
 208                ifld(n)=ifld(n)+gref(j)
 209                n=n+1
 210              enddo
 211            else
 212              ifld(n:n+glen(j)-1)=gref(j)
 213              n=n+glen(j)
 214            endif
 215            iofst=iofst+(gwidth(j)*glen(j))
 216          enddo
 217       elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then
 218          ! missing values included
 219          allocate(ifldmiss(ndpts))
 220          !ifldmiss=0
 221          n=1
 222          non=1
 223          do j=1,ngroups
 224            !print *,'SAGNGP ',j,gwidth(j),glen(j),gref(j)
 225            if (gwidth(j).ne.0) then
 226              msng1=(2**gwidth(j))-1
 227              msng2=msng1-1
 228              call gbytes(cpack,ifld(n),iofst,gwidth(j),0,glen(j))
 229              iofst=iofst+(gwidth(j)*glen(j))
 230              do k=1,glen(j)
 231                if (ifld(n).eq.msng1) then
 232                   ifldmiss(n)=1
 233                elseif (idrstmpl(7).eq.2.AND.ifld(n).eq.msng2) then
 234                   ifldmiss(n)=2
 235                else
 236                   ifldmiss(n)=0
 237                   ifld(non)=ifld(n)+gref(j)
 238                   non=non+1
 239                endif
 240                n=n+1
 241              enddo
 242            else
 243              msng1=(2**nbitsgref)-1
 244              msng2=msng1-1
 245              if (gref(j).eq.msng1) then
 246                 ifldmiss(n:n+glen(j)-1)=1
 247                 !ifld(n:n+glen(j)-1)=0
 248              elseif (idrstmpl(7).eq.2.AND.gref(j).eq.msng2) then
 249                 ifldmiss(n:n+glen(j)-1)=2
 250                 !ifld(n:n+glen(j)-1)=0
 251              else
 252                 ifldmiss(n:n+glen(j)-1)=0
 253                 ifld(non:non+glen(j)-1)=gref(j)
 254                 non=non+glen(j)
 255              endif
 256              n=n+glen(j)
 257            endif
 258          enddo
 259       endif
 260       !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts)
 261
 262       if ( allocated(gref) ) deallocate(gref)
 263       if ( allocated(gwidth) ) deallocate(gwidth)
 264       if ( allocated(glen) ) deallocate(glen)
 265 !
 266 !  If using spatial differences, add overall min value, and
 267 !  sum up recursively
 268 !
 269       if (idrsnum.eq.3) then         ! spatial differencing
 270          if (idrstmpl(17).eq.1) then      ! first order
 271             ifld(1)=ival1
 272             if ( idrstmpl(7).eq.0 ) then        ! no missing values
 273                itemp=ndpts
 274             else
 275                itemp=non-1
 276             endif
 277             do n=2,itemp
 278                ifld(n)=ifld(n)+minsd
 279                ifld(n)=ifld(n)+ifld(n-1)
 280             enddo
 281          elseif (idrstmpl(17).eq.2) then    ! second order
 282             ifld(1)=ival1
 283             ifld(2)=ival2
 284             if ( idrstmpl(7).eq.0 ) then        ! no missing values
 285                itemp=ndpts
 286             else
 287                itemp=non-1
 288             endif
 289             do n=3,itemp
 290                ifld(n)=ifld(n)+minsd
 291                ifld(n)=ifld(n)+(2*ifld(n-1))-ifld(n-2)
 292             enddo
 293          endif
 294       !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts)
 295       endif
 296 !
 297 !  Scale data back to original form
 298 !
 299       !print *,'SAGT: ',ref,bscale,dscale
 300       if ( idrstmpl(7).eq.0 ) then        ! no missing values
 301          do n=1,ndpts
 302             fld(n)=((real(ifld(n))*bscale)+ref)*dscale
 303             !write(78,*)'SAG ',n,fld(n),ifld(n),bscale,ref,1./dscale
 304          enddo
 305       elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then
 306          ! missing values included
 307          non=1
 308          do n=1,ndpts
 309             if ( ifldmiss(n).eq.0 ) then
 310                fld(n)=((real(ifld(non))*bscale)+ref)*dscale
 311                !print *,'SAG ',n,fld(n),ifld(non),bscale,ref,dscale
 312                non=non+1
 313             elseif ( ifldmiss(n).eq.1 ) then
 314                fld(n)=rmiss1
 315             elseif ( ifldmiss(n).eq.2 ) then
 316                fld(n)=rmiss2
 317             endif
 318          enddo
 319          if ( allocated(ifldmiss) ) deallocate(ifldmiss)
 320       endif
 321
 322       if ( allocated(ifld) ) deallocate(ifld)
 323
 324       return
 325       end