source: lmdz_wrf/trunk/WRFV3/external/io_grib2/g2lib/comunpack.F

Last change on this file was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 10.5 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.