1 subroutine gf_getfld(cgrib,lcgrib,ifldnum,unpack,expand,gfld,ierr)
2 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
4 ! SUBPROGRAM: gf_getfld
5 ! PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-05-26
7 ! ABSTRACT: This subroutine returns the Grid Definition, Product Definition,
8 ! Bit-map ( if applicable ), and the unpacked data for a given data
9 ! field. All of the information returned is stored in a derived
10 ! type variable, gfld. Gfld is of type gribfield, which is defined
11 ! in module grib_mod, so users of this routine will need to include
12 ! the line "USE GRIB_MOD" in their calling routine. Each component of the
13 ! gribfield type is described in the OUTPUT ARGUMENT LIST section below.
15 ! Since there can be multiple data fields packed into a GRIB2
16 ! message, the calling routine indicates which field is being requested
17 ! with the ifldnum argument.
19 ! PROGRAM HISTORY LOG:
21 ! 2002-01-24 Gilbert - Changed to pass back derived type gribfield
22 ! variable through argument list, instead of
23 ! having many different arguments.
24 ! 2004-05-20 Gilbert - Added check to see if previous a bit-map is specified,
27 ! USAGE: CALL gf_getfld(cgrib,lcgrib,ifldnum,unpack,expand,gfld,ierr)
28 ! INPUT ARGUMENT LIST:
29 ! cgrib - Character array that contains the GRIB2 message
30 ! lcgrib - Length (in bytes) of GRIB message array cgrib.
31 ! ifldnum - Specifies which field in the GRIB2 message to return.
32 ! unpack - Logical value indicating whether to unpack bitmap/data
33 ! .true. = unpack bitmap and data values
34 ! .false. = do not unpack bitmap and data values
35 ! expand - Boolean value indicating whether the data points should be
36 ! expanded to the correspond grid, if a bit-map is present.
37 ! 1 = if possible, expand data field to grid, inserting zero
38 ! values at gridpoints that are bitmapped out.
40 ! 0 = do not expand data field, leaving it an array of
41 ! consecutive data points for each "1" in the bitmap.
42 ! This argument is ignored if unpack == 0 OR if the
43 ! returned field does not contain a bit-map.
45 ! OUTPUT ARGUMENT LIST:
46 ! gfld - derived type gribfield ( defined in module grib_mod )
47 ! ( NOTE: See Remarks Section )
48 ! gfld%version = GRIB edition number ( currently 2 )
49 ! gfld%discipline = Message Discipline ( see Code Table 0.0 )
50 ! gfld%idsect() = Contains the entries in the Identification
51 ! Section ( Section 1 )
52 ! This element is actually a pointer to an array
53 ! that holds the data.
54 ! gfld%idsect(1) = Identification of originating Centre
55 ! ( see Common Code Table C-1 )
56 ! 7 - US National Weather Service
57 ! gfld%idsect(2) = Identification of originating Sub-centre
58 ! gfld%idsect(3) = GRIB Master Tables Version Number
59 ! ( see Code Table 1.0 )
61 ! 1 - Initial operational version number
62 ! gfld%idsect(4) = GRIB Local Tables Version Number
63 ! ( see Code Table 1.1 )
64 ! 0 - Local tables not used
65 ! 1-254 - Number of local tables version used
66 ! gfld%idsect(5) = Significance of Reference Time (Code Table 1.2)
68 ! 1 - Start of forecast
69 ! 2 - Verifying time of forecast
70 ! 3 - Observation time
71 ! gfld%idsect(6) = Year ( 4 digits )
72 ! gfld%idsect(7) = Month
73 ! gfld%idsect(8) = Day
74 ! gfld%idsect(9) = Hour
75 ! gfld%idsect(10) = Minute
76 ! gfld%idsect(11) = Second
77 ! gfld%idsect(12) = Production status of processed data
78 ! ( see Code Table 1.3 )
79 ! 0 - Operational products
80 ! 1 - Operational test products
81 ! 2 - Research products
82 ! 3 - Re-analysis products
83 ! gfld%idsect(13) = Type of processed data ( see Code Table 1.4 )
84 ! 0 - Analysis products
85 ! 1 - Forecast products
86 ! 2 - Analysis and forecast products
87 ! 3 - Control forecast products
88 ! 4 - Perturbed forecast products
89 ! 5 - Control and perturbed forecast products
90 ! 6 - Processed satellite observations
91 ! 7 - Processed radar observations
92 ! gfld%idsectlen = Number of elements in gfld%idsect().
93 ! gfld%local() = Pointer to character array containing contents
94 ! of Local Section 2, if included
95 ! gfld%locallen = length of array gfld%local()
96 ! gfld%ifldnum = field number within GRIB message
97 ! gfld%griddef = Source of grid definition (see Code Table 3.0)
98 ! 0 - Specified in Code table 3.1
99 ! 1 - Predetermined grid Defined by originating centre
100 ! gfld%ngrdpts = Number of grid points in the defined grid.
101 ! gfld%numoct_opt = Number of octets needed for each
102 ! additional grid points definition.
103 ! Used to define number of
104 ! points in each row ( or column ) for
106 ! = 0, if using regular grid.
107 ! gfld%interp_opt = Interpretation of list for optional points
108 ! definition. (Code Table 3.11)
109 ! gfld%igdtnum = Grid Definition Template Number (Code Table 3.1)
110 ! gfld%igdtmpl() = Contains the data values for the specified Grid
111 ! Definition Template ( NN=gfld%igdtnum ). Each
112 ! element of this integer array contains an entry (in
113 ! the order specified) of Grid Defintion Template 3.NN
114 ! This element is actually a pointer to an array
115 ! that holds the data.
116 ! gfld%igdtlen = Number of elements in gfld%igdtmpl(). i.e. number of
117 ! entries in Grid Defintion Template 3.NN
118 ! ( NN=gfld%igdtnum ).
119 ! gfld%list_opt() = (Used if gfld%numoct_opt .ne. 0) This array
120 ! contains the number of grid points contained in
121 ! each row ( or column ). (part of Section 3)
122 ! This element is actually a pointer to an array
123 ! that holds the data. This pointer is nullified
124 ! if gfld%numoct_opt=0.
125 ! gfld%num_opt = (Used if gfld%numoct_opt .ne. 0) The number of entries
126 ! in array ideflist. i.e. number of rows ( or columns )
127 ! for which optional grid points are defined. This value
128 ! is set to zero, if gfld%numoct_opt=0.
129 ! gfdl%ipdtnum = Product Definition Template Number (see Code Table 4.0)
130 ! gfld%ipdtmpl() = Contains the data values for the specified Product
131 ! Definition Template ( N=gfdl%ipdtnum ). Each element
132 ! of this integer array contains an entry (in the
133 ! order specified) of Product Defintion Template 4.N.
134 ! This element is actually a pointer to an array
135 ! that holds the data.
136 ! gfld%ipdtlen = Number of elements in gfld%ipdtmpl(). i.e. number of
137 ! entries in Product Defintion Template 4.N
138 ! ( N=gfdl%ipdtnum ).
139 ! gfld%coord_list() = Real array containing floating point values
140 ! intended to document the vertical discretisation
141 ! associated to model data on hybrid coordinate
142 ! vertical levels. (part of Section 4)
143 ! This element is actually a pointer to an array
144 ! that holds the data.
145 ! gfld%num_coord = number of values in array gfld%coord_list().
146 ! gfld%ndpts = Number of data points unpacked and returned.
147 ! gfld%idrtnum = Data Representation Template Number
148 ! ( see Code Table 5.0)
149 ! gfld%idrtmpl() = Contains the data values for the specified Data
150 ! Representation Template ( N=gfld%idrtnum ). Each
151 ! element of this integer array contains an entry
152 ! (in the order specified) of Product Defintion
154 ! This element is actually a pointer to an array
155 ! that holds the data.
156 ! gfld%idrtlen = Number of elements in gfld%idrtmpl(). i.e. number
157 ! of entries in Data Representation Template 5.N
158 ! ( N=gfld%idrtnum ).
159 ! gfld%unpacked = logical value indicating whether the bitmap and
160 ! data values were unpacked. If false,
161 ! gfld%bmap and gfld%fld pointers are nullified.
162 ! gfld%expanded = Logical value indicating whether the data field
163 ! was expanded to the grid in the case where a
164 ! bit-map is present. If true, the data points in
165 ! gfld%fld match the grid points and zeros were
166 ! inserted at grid points where data was bit-mapped
167 ! out. If false, the data values in gfld%fld were
168 ! not expanded to the grid and are just a consecutive
169 ! array of data points corresponding to each value of
171 ! gfld%ibmap = Bitmap indicator ( see Code Table 6.0 )
172 ! 0 = bitmap applies and is included in Section 6.
173 ! 1-253 = Predefined bitmap applies
174 ! 254 = Previously defined bitmap applies to this field
175 ! 255 = Bit map does not apply to this product.
176 ! gfld%bmap() = Logical*1 array containing decoded bitmap,
177 ! if ibmap=0 or ibap=254. Otherwise nullified.
178 ! This element is actually a pointer to an array
179 ! that holds the data.
180 ! gfld%fld() = Array of gfld%ndpts unpacked data points.
181 ! This element is actually a pointer to an array
182 ! that holds the data.
183 ! ierr - Error return code.
185 ! 1 = Beginning characters "GRIB" not found.
186 ! 2 = GRIB message is not Edition 2.
187 ! 3 = The data field request number was not positive.
188 ! 4 = End string "7777" found, but not where expected.
189 ! 6 = GRIB message did not contain the requested number of
191 ! 7 = End string "7777" not found at end of message.
192 ! 8 = Unrecognized Section encountered.
193 ! 9 = Data Representation Template 5.NN not yet implemented.
194 ! 15 = Error unpacking Section 1.
195 ! 16 = Error unpacking Section 2.
196 ! 10 = Error unpacking Section 3.
197 ! 11 = Error unpacking Section 4.
198 ! 12 = Error unpacking Section 5.
199 ! 13 = Error unpacking Section 6.
200 ! 14 = Error unpacking Section 7.
201 ! 17 = Previous bitmap specified, but none exists.
203 ! REMARKS: Note that derived type gribfield contains pointers to many
204 ! arrays of data. The memory for these arrays is allocated
205 ! when the values in the arrays are set, to help minimize
206 ! problems with array overloading. Because of this users
207 ! are encouraged to free up this memory, when it is no longer
208 ! needed, by an explicit call to subroutine gf_free.
209 ! ( i.e. CALL GF_FREE(GFLD) )
211 ! Subroutine gb_info can be used to first determine
212 ! how many data fields exist in a given GRIB message.
214 ! REMARKS2: It may not always be possible to expand a bit-mapped data field.
215 ! If a pre-defined bit-map is used and not included in the GRIB2
216 ! message itself, this routine would not have the necessary
217 ! information to expand the data. In this case, gfld%expanded would
218 ! would be set to 0 (false), regardless of the value of input
222 ! LANGUAGE: Fortran 90
228 character(len=1),intent(in) :: cgrib(lcgrib)
229 integer,intent(in) :: lcgrib,ifldnum
230 logical,intent(in) :: unpack,expand
231 type(gribfield),intent(out) :: gfld
232 integer,intent(out) :: ierr
233 ! integer,intent(out) :: igds(*),igdstmpl(*),ideflist(*)
234 ! integer,intent(out) :: ipdsnum,ipdstmpl(*)
235 ! integer,intent(out) :: idrsnum,idrstmpl(*)
236 ! integer,intent(out) :: ndpts,ibmap,idefnum,numcoord
237 ! logical*1,intent(out) :: bmap(*)
238 ! real,intent(out) :: fld(*),coordlist(*)
240 character(len=4),parameter :: grib='GRIB',c7777='7777'
241 character(len=4) :: ctemp
242 real,pointer,dimension(:) :: newfld
243 integer:: listsec0(2),igds(5)
244 integer iofst,ibeg,istart
246 logical*1,pointer,dimension(:) :: bmpsave
247 logical have3,have4,have5,have6,have7
250 subroutine gf_unpack1(cgrib,lcgrib,iofst,ids,idslen,ierr)
251 character(len=1),intent(in) :: cgrib(lcgrib)
252 integer,intent(in) :: lcgrib
253 integer,intent(inout) :: iofst
254 integer,pointer,dimension(:) :: ids
255 integer,intent(out) :: ierr,idslen
256 end subroutine gf_unpack1
257 subroutine gf_unpack2(cgrib,lcgrib,iofst,lencsec2,csec2,ierr)
258 character(len=1),intent(in) :: cgrib(lcgrib)
259 integer,intent(in) :: lcgrib
260 integer,intent(inout) :: iofst
261 integer,intent(out) :: lencsec2
262 integer,intent(out) :: ierr
263 character(len=1),pointer,dimension(:) :: csec2
264 end subroutine gf_unpack2
265 subroutine gf_unpack3(cgrib,lcgrib,iofst,igds,igdstmpl,
266 & mapgridlen,ideflist,idefnum,ierr)
267 character(len=1),intent(in) :: cgrib(lcgrib)
268 integer,intent(in) :: lcgrib
269 integer,intent(inout) :: iofst
270 integer,pointer,dimension(:) :: igdstmpl,ideflist
271 integer,intent(out) :: igds(5)
272 integer,intent(out) :: ierr,idefnum
273 end subroutine gf_unpack3
274 subroutine gf_unpack4(cgrib,lcgrib,iofst,ipdsnum,ipdstmpl,
275 & mappdslen,coordlist,numcoord,ierr)
276 character(len=1),intent(in) :: cgrib(lcgrib)
277 integer,intent(in) :: lcgrib
278 integer,intent(inout) :: iofst
279 real,pointer,dimension(:) :: coordlist
280 integer,pointer,dimension(:) :: ipdstmpl
281 integer,intent(out) :: ipdsnum
282 integer,intent(out) :: ierr,numcoord
283 end subroutine gf_unpack4
284 subroutine gf_unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,
285 & idrstmpl,mapdrslen,ierr)
286 character(len=1),intent(in) :: cgrib(lcgrib)
287 integer,intent(in) :: lcgrib
288 integer,intent(inout) :: iofst
289 integer,intent(out) :: ndpts,idrsnum
290 integer,pointer,dimension(:) :: idrstmpl
291 integer,intent(out) :: ierr
292 end subroutine gf_unpack5
293 subroutine gf_unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
294 character(len=1),intent(in) :: cgrib(lcgrib)
295 integer,intent(in) :: lcgrib,ngpts
296 integer,intent(inout) :: iofst
297 integer,intent(out) :: ibmap
298 integer,intent(out) :: ierr
299 logical*1,pointer,dimension(:) :: bmap
300 end subroutine gf_unpack6
301 subroutine gf_unpack7(cgrib,lcgrib,iofst,igdsnum,igdstmpl,
302 & idrsnum,idrstmpl,ndpts,fld,ierr)
303 character(len=1),intent(in) :: cgrib(lcgrib)
304 integer,intent(in) :: lcgrib,ndpts,idrsnum,igdsnum
305 integer,intent(inout) :: iofst
306 integer,pointer,dimension(:) :: idrstmpl,igdstmpl
307 integer,intent(out) :: ierr
308 real,pointer,dimension(:) :: fld
309 end subroutine gf_unpack7
320 nullify(gfld%list_opt,gfld%igdtmpl,gfld%ipdtmpl)
321 nullify(gfld%coord_list,gfld%idrtmpl,gfld%bmap,gfld%fld)
323 ! Check for valid request number
325 if (ifldnum.le.0) then
326 print *,'gf_getfld: Request for field number must be positive.'
331 ! Check for beginning of GRIB message in the first 100 bytes
335 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
336 if (ctemp.eq.grib ) then
341 if (istart.eq.0) then
342 print *,'gf_getfld: Beginning characters GRIB not found.'
347 ! Unpack Section 0 - Indicator Section
350 call g2lib_gbyte(cgrib,listsec0(1),iofst,8) ! Discipline
352 call g2lib_gbyte(cgrib,listsec0(2),iofst,8) ! GRIB edition number
355 call g2lib_gbyte(cgrib,lengrib,iofst,32) ! Length of GRIB message
360 ! Currently handles only GRIB Edition 2.
362 if (listsec0(2).ne.2) then
363 print *,'gf_getfld: can only decode GRIB edition 2.'
368 ! Loop through the remaining sections keeping track of the
369 ! length of each. Also keep the latest Grid Definition Section info.
370 ! Unpack the requested field number.
373 ! Check to see if we are at end of GRIB message
374 ctemp=cgrib(ipos)//cgrib(ipos+1)//cgrib(ipos+2)//cgrib(ipos+3)
375 if (ctemp.eq.c7777 ) then
377 ! If end of GRIB message not where expected, issue error
378 if (ipos.ne.(istart+lengrib)) then
379 print *,'gf_getfld: "7777" found, but not where expected.'
385 ! Get length of Section and Section number
387 call g2lib_gbyte(cgrib,lensec,iofst,32) ! Get Length of Section
389 call g2lib_gbyte(cgrib,isecnum,iofst,8) ! Get Section number
391 !print *,' lensec= ',lensec,' secnum= ',isecnum
393 ! Check to see if section number is valid
395 if ( (isecnum.lt.1).OR.(isecnum.gt.7) ) then
396 print *,'gf_getfld: Unrecognized Section Encountered=',isecnum
401 ! If found Section 1, decode elements in Identification Section
403 if (isecnum.eq.1) then
404 iofst=iofst-40 ! reset offset to beginning of section
405 call gf_unpack1(cgrib,lcgrib,iofst,gfld%idsect,
406 & gfld%idsectlen,jerr)
413 ! If found Section 2, Grab local section
414 ! Save in case this is the latest one before the requested field.
416 if (isecnum.eq.2) then
417 iofst=iofst-40 ! reset offset to beginning of section
418 if (associated(gfld%local)) deallocate(gfld%local)
419 call gf_unpack2(cgrib,lcgrib,iofst,gfld%locallen,
427 ! If found Section 3, unpack the GDS info using the
428 ! appropriate template. Save in case this is the latest
429 ! grid before the requested field.
431 if (isecnum.eq.3) then
432 iofst=iofst-40 ! reset offset to beginning of section
433 if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl)
434 if (associated(gfld%list_opt)) deallocate(gfld%list_opt)
435 call gf_unpack3(cgrib,lcgrib,iofst,igds,gfld%igdtmpl,
436 & gfld%igdtlen,gfld%list_opt,gfld%num_opt,jerr)
441 gfld%numoct_opt=igds(3)
442 gfld%interp_opt=igds(4)
450 ! If found Section 4, check to see if this field is the
453 if (isecnum.eq.4) then
455 if (numfld.eq.ifldnum) then
456 gfld%discipline=listsec0(1)
457 gfld%version=listsec0(2)
460 gfld%expanded=.false.
461 iofst=iofst-40 ! reset offset to beginning of section
462 call gf_unpack4(cgrib,lcgrib,iofst,gfld%ipdtnum,
463 & gfld%ipdtmpl,gfld%ipdtlen,gfld%coord_list,
464 & gfld%num_coord,jerr)
474 ! If found Section 5, check to see if this field is the
477 if ((isecnum.eq.5).and.(numfld.eq.ifldnum)) then
478 iofst=iofst-40 ! reset offset to beginning of section
479 call gf_unpack5(cgrib,lcgrib,iofst,gfld%ndpts,gfld%idrtnum,
480 & gfld%idrtmpl,gfld%idrtlen,jerr)
489 ! If found Section 6, Unpack bitmap.
490 ! Save in case this is the latest
491 ! bitmap before the requested field.
493 if (isecnum.eq.6) then
494 if (unpack) then ! unpack bitmap
495 iofst=iofst-40 ! reset offset to beginning of section
496 bmpsave=>gfld%bmap ! save pointer to previous bitmap
497 call gf_unpack6(cgrib,lcgrib,iofst,gfld%ngrdpts,gfld%ibmap,
501 if (gfld%ibmap .eq. 254) then ! use previously specified bitmap
502 if ( associated(bmpsave) ) then
505 print *,'gf_getfld: Previous bit-map specified,',
506 & ' but none exists,'
511 if ( associated(bmpsave) ) deallocate(bmpsave)
517 else ! do not unpack bitmap
518 call g2lib_gbyte(cgrib,gfld%ibmap,iofst,8) ! Get BitMap Indicator
523 ! If found Section 7, check to see if this field is the
526 if ((isecnum.eq.7).and.(numfld.eq.ifldnum).and.unpack) then
527 iofst=iofst-40 ! reset offset to beginning of section
528 call gf_unpack7(cgrib,lcgrib,iofst,gfld%igdtnum,
529 & gfld%igdtmpl,gfld%idrtnum,
530 & gfld%idrtmpl,gfld%ndpts,
534 ! If bitmap is used with this field, expand data field
535 ! to grid, if possible.
536 if ( gfld%ibmap .ne. 255 .AND. associated(gfld%bmap) ) then
538 allocate(newfld(gfld%ngrdpts))
539 !newfld(1:gfld%ngrdpts)=0.0
540 !newfld=unpack(gfld%fld,gfld%bmap,newfld)
543 if ( gfld%bmap(j) ) then
544 newfld(j)=gfld%fld(n)
550 deallocate(gfld%fld);
554 gfld%expanded=.false.
560 print *,'gf_getfld: return from gf_unpack7 = ',jerr
566 ! Check to see if we read pass the end of the GRIB
567 ! message and missed the terminator string '7777'.
569 ipos=ipos+lensec ! Update beginning of section pointer
570 if (ipos.gt.(istart+lengrib)) then
571 print *,'gf_getfld: "7777" not found at end of GRIB message.'
576 ! If unpacking requested, return when all sections have been
579 if (unpack.and.have3.and.have4.and.have5.and.have6.and.have7)
582 ! If unpacking is not requested, return when sections
583 ! 3 through 6 have been processed
585 if ((.NOT.unpack).and.have3.and.have4.and.have5.and.have6)
591 ! If exited from above loop, the end of the GRIB message was reached
592 ! before the requested field was found.
594 print *,'gf_getfld: GRIB message contained ',numlocal,
595 & ' different fields.'
596 print *,'gf_getfld: The request was for the ',ifldnum,