Context Navigation

comunpack.F

Last change on this file was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
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.

Download in other formats:

Original Format