| 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 |
|---|
