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