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

Rev	Line
[2759]	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