emcsfc_snow2mdl 1.14.0
Loading...
Searching...
No Matches
grib_utils.F90
Go to the documentation of this file.
1
4
23 subroutine grib_check(file_name, isgrib)
24
25 implicit none
26
27 character*(*), intent(in) :: file_name
28 integer, parameter :: iunit=11
29 integer :: istat, iseek, mseek, lskip, lgrib, version
30 integer, intent(out) :: isgrib
31
32 print*,"- CHECK FILE TYPE OF: ", trim(file_name)
33 call baopenr (iunit, file_name, istat)
34
35 if (istat /= 0) then
36 print*,'- FATAL ERROR: BAD FILE OPEN. ISTAT IS ',istat
37 call w3tage('SNOW2MDL')
38 call errexit(40)
39 end if
40
41 iseek = 0
42 mseek = 64
43 call skgb2(iunit, iseek, mseek, lskip, lgrib, version)
44
45 call baclose(iunit, istat)
46
47 if (lgrib > 0) then
48 isgrib = version
49 if (isgrib == 1) print*,"- FILE IS GRIB1"
50 if (isgrib == 2) print*,"- FILE IS GRIB2"
51 else
52 isgrib = 0
53 print*,"- FILE IS BINARY"
54 endif
55
56 return
57
58 end subroutine grib_check
59
75 SUBROUTINE skgb2(LUGB,ISEEK,MSEEK,LSKIP,LGRIB,I1)
76 implicit none
77 INTEGER, INTENT( IN) :: LUGB, ISEEK, MSEEK
78 INTEGER, INTENT(OUT) :: LSKIP, LGRIB, I1
79 INTEGER, PARAMETER :: LSEEK=128
80 INTEGER :: K, KZ, KS, KG, KN, KM, I4, K4
81 CHARACTER Z(LSEEK)
82 CHARACTER Z4(4)
83! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84 i1=0
85 lgrib=0
86 ks=iseek
87 kn=min(lseek,mseek)
88 kz=lseek
89! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90! LOOP UNTIL GRIB MESSAGE IS FOUND
91 DO WHILE(lgrib.EQ.0.AND.kn.GE.8.AND.kz.EQ.lseek)
92! READ PARTIAL SECTION
93 CALL baread(lugb,ks,kn,kz,z)
94 km=kz-8+1
95 k=0
96! LOOK FOR 'GRIB...1' IN PARTIAL SECTION
97 DO WHILE(lgrib.EQ.0.AND.k.LT.km)
98 CALL gbytec(z,i4,(k+0)*8,4*8)
99 CALL gbytec(z,i1,(k+7)*8,1*8)
100 IF(i4.EQ.1196575042.AND.(i1.EQ.1.OR.i1.EQ.2)) THEN
101! LOOK FOR '7777' AT END OF GRIB MESSAGE
102 IF (i1.EQ.1) CALL gbytec(z,kg,(k+4)*8,3*8)
103 IF (i1.EQ.2) CALL gbytec(z,kg,(k+12)*8,4*8)
104 CALL baread(lugb,ks+k+kg-4,4,k4,z4)
105 IF(k4.EQ.4) THEN
106 CALL gbytec(z4,i4,0,4*8)
107 IF(i4.EQ.926365495) THEN
108! GRIB MESSAGE FOUND
109 lskip=ks+k
110 lgrib=kg
111 ENDIF
112 ENDIF
113 ENDIF
114 k=k+1
115 ENDDO
116 ks=ks+km
117 kn=min(lseek,iseek+mseek-ks)
118 ENDDO
119! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120 RETURN
121 END subroutine skgb2
122
139 subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
140
141 implicit none
142
143 integer, intent(in ) :: igdtnum, igdtlen, igdstmpl(igdtlen)
144 integer, intent( out) :: kgds(200), ni, nj
145 integer :: iscale
146
147 real, intent( out) :: res
148
149 kgds=0
150
151 if (igdtnum.eq.0) then ! lat/lon grid
152
153 iscale=igdstmpl(10)*igdstmpl(11)
154 if (iscale == 0) iscale = 1e6
155 kgds(1)=0 ! oct 6
156 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
157 ni = kgds(2)
158 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
159 nj = kgds(3)
160 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
161 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
162
163 kgds(6)=0 ! oct 17, resolution and component flags
164 if (igdstmpl(1)==2 ) kgds(6)=64
165 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
166 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
167
168 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point
169 kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point
170 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, di
171 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj
172
173 kgds(11) = 0 ! oct 28, scan mode
174 if (btest(igdstmpl(19),7)) kgds(11) = 128
175 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
176 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
177
178 kgds(12)=0 ! octs 29-32, reserved
179 kgds(19)=0 ! oct 4, # vert coordinate parameters
180 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
181
182 res = float(kgds(9)) / 1000.0 * 111.0
183
184 elseif (igdtnum.eq.40) then ! Gaussian Lat/Lon grid
185
186 iscale=igdstmpl(10)*igdstmpl(11)
187 if (iscale==0) iscale=1e6
188 kgds(1)=4 ! oct 6
189 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
190 ni = kgds(2)
191 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
192 nj = kgds(3)
193 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
194 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
195
196 kgds(6)=0 ! oct 17, resolution and component flags
197 if (igdstmpl(1)==2 ) kgds(6)=64
198 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
199 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
200
201 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point
202 kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point
203 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, Di
204 kgds(10)=igdstmpl(18) ! octs 26-27, Number of parallels
205
206 kgds(11) = 0 ! oct 28, scan mode
207 if (btest(igdstmpl(19),7)) kgds(11) = 128
208 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
209 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
210
211 kgds(12)=0 ! octs 29-32, reserved
212 kgds(19)=0 ! oct 4, # vert coordinate parameters
213 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
214
215 res = float(kgds(9)) / 1000.0 * 111.0
216
217 elseif (igdtnum.eq.20) then ! Polar Stereographic Grid
218
219 iscale=1e6
220 kgds(1)=5 ! oct 6, data representation type, polar
221 kgds(2)=igdstmpl(8) ! octs 7-8, nx
222 ni = kgds(2)
223 kgds(3)=igdstmpl(9) ! octs 8-10, ny
224 nj = kgds(3)
225 kgds(4)=nint(float(igdstmpl(10))/float(iscale)*1000.) ! octs 11-13, lat of 1st grid point
226 kgds(5)=nint(float(igdstmpl(11))/float(iscale)*1000.) ! octs 14-16, lon of 1st grid point
227
228 kgds(6)=0 ! oct 17, resolution and component flags
229 if (igdstmpl(1) >= 2 .or. igdstmpl(1) <= 5) kgds(6)=64
230 if (igdstmpl(1) == 7) kgds(6)=64
231 if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128
232 if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
233
234 kgds(7)=nint(float(igdstmpl(14))/float(iscale)*1000.) ! octs 18-20, lon of orientation
235 kgds(8)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 21-23, dx
236 kgds(9)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 24-26, dy
237
238 kgds(10)=0 ! oct 27, projection center flag
239 if (btest(igdstmpl(17),1)) kgds(10) = 128
240
241 kgds(11) = 0 ! oct 28, scan mode
242 if (btest(igdstmpl(18),7)) kgds(11) = 128
243 if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) + 64
244 if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) + 32
245
246 kgds(19)=0 ! oct 4, # vert coordinate parameters
247 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
248
249 res = 0.5 * float(kgds(8)+kgds(9)) / 1000.
250
251 elseif (igdtnum.eq.1) then ! Rotated Lat/Lon grid
252
253 if (btest(igdstmpl(19),2)) then ! e-stagger, bit 6 of scan mode is '1'
254
255 iscale=igdstmpl(10)*igdstmpl(11)
256 if (iscale == 0) iscale = 1e6
257 kgds(1)=203 ! oct 6, "E" grid
258 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
259 ni = kgds(2)
260 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
261 nj = kgds(3)
262 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
263 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
264
265 kgds(6)=0 ! oct 17, resolution and component flags
266 if (igdstmpl(1)==2 ) kgds(6)=64
267 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
268 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
269
270 kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000 ! octs 18-20, Lat of cent of rotation
271 kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23, Lon of cent of rotation
272 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*500.) ! octs 24-25, Di
273 ! Note!! grib 2 convention twice grib 1
274 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, Dj
275
276 kgds(11) = 0 ! oct 28, scan mode
277 if (btest(igdstmpl(19),7)) kgds(11) = 128
278 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
279 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
280
281 kgds(12)=0 ! octs 29-32, reserved
282 kgds(19)=0 ! oct 4, # vert coordinate parameters
283 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
284
285 res = sqrt( (float(kgds(9)) / 1000.0)**2 + &
286 (float(kgds(10)) / 1000.0)**2 )
287 res = res * 111.0
288
289 else ! b-stagger
290
291 iscale=igdstmpl(10)*igdstmpl(11)
292 if (iscale == 0) iscale = 1e6
293 kgds(1)=205 ! oct 6, rotated lat/lon for Non-E Stagger grid
294 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
295 ni = kgds(2)
296 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
297 nj = kgds(3)
298 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
299 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
300
301 kgds(6)=0 ! oct 17, resolution and component flags
302 if (igdstmpl(1)==2 ) kgds(6)=64
303 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
304 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
305
306 kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000 ! octs 18-20, Lat of cent of rotation
307 kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23, Lon of cent of rotation
308 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, Di
309 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, Dj
310
311 kgds(11) = 0 ! oct 28, scan mode
312 if (btest(igdstmpl(19),7)) kgds(11) = 128
313 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
314 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
315
316 kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of last grid point
317 kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of last grid point
318
319 kgds(19)=0 ! oct 4, # vert coordinate parameters
320 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
321
322 res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
323 * 0.5 * 111.0
324
325 endif
326
327 else
328
329 print*,'- FATAL ERROR CONVERTING TO GRIB2 GDT.'
330 print*,'- UNRECOGNIZED GRID TYPE.'
331 call w3tage('SNOW2MDL')
332 call errexit(50)
333
334 endif
335
336 end subroutine gdt_to_gds
337
350 subroutine grib2_check (kgds, igdstmplen)
351 implicit none
352
353 integer, intent(in) :: kgds(200)
354 integer, intent( out) :: igdstmplen
355
356 select case (kgds(1))
357 case(4) ! gaussian
358 igdstmplen = 19
359 case(203, 205) ! rotated lat/lon "B" or "E" stagger
360 igdstmplen = 22
361 case default
362 print*,'- FATAL ERROR IN ROUTINE GRIB2_CHECK.'
363 print*,'- UNRECOGNIZED GRID TYPE.'
364 call w3tage('SNOW2MDL')
365 call errexit(47)
366 end select
367
368 end subroutine grib2_check
369
397 subroutine init_grib2(century, year, month, day, hour, kgds, &
398 lat11, latlast, lon11, lonlast, &
399 listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl, &
400 igdstmplen, idefnum, ideflist, ngrdpts)
401 implicit none
402
403 integer, intent(in ) :: century, year, month, day, hour
404 integer, intent(in ) :: kgds(200), igdstmplen
405 integer, intent( out) :: igds(5)
406 integer, intent( out) :: listsec0(2)
407 integer, intent( out) :: listsec1(13)
408 integer, intent( out) :: ipdstmpl(15), ipdsnum
409 integer, intent( out) :: igdstmpl(igdstmplen)
410 integer, intent( out) :: idefnum, ideflist
411 integer, intent( out) :: ngrdpts
412
413 real, intent(in ) :: lat11, latlast, lon11, lonlast
414 real :: scale
415
416! Section 0
417
418 listsec0(1)=0 ! discipline, meteorological fields
419 listsec0(2)=2 ! grib version 2
420
421! Section 1
422
423 listsec1(1)=7 ! id of center (ncep)
424 listsec1(2)=4 ! subcenter (emc)
425 listsec1(3)=8 ! master table version number. wgrib2 does not recognize later tables
426 listsec1(4)= 0 ! local table not used
427 listsec1(5)= 0 ! signif of ref time - analysis
428 if (year == 100) then
429 listsec1(6)=century*100 + year
430 else
431 listsec1(6)=(century-1)*100 + year
432 endif
433 listsec1(7)=month
434 listsec1(8)=day
435 listsec1(9)=hour
436 listsec1(10:11)=0 ! minutes/secs
437 listsec1(12)=0 ! production status of data - ops products
438 listsec1(13)=0 ! type of processed products - analysis
439
440! Section 2 - not used
441
442! Section 3 - grid description section
443
444 if (kgds(1) == 4) then ! gaussian
445
446 igdstmpl(1)=5 ! oct 15; shape of the earth, wgs84
447 igdstmpl(2)=255 ! oct 16; scale factor of radius of spherical earth, not used.
448 igdstmpl(3)=-1 ! octs 17-20; scale value of radius of spherical earth, not used.
449 igdstmpl(4)=255 ! oct 21; scale factor of major axis of elliptical earth, not used.
450 igdstmpl(5)=-1 ! octs 22-25; scaled value of major axis of elliptical earth, not used.
451 igdstmpl(6)=255 ! oct 26; scale factor of minor axis of elliptical earth, not used.
452 igdstmpl(7)=-1 ! octs 27-30; scaled value of minor axis of elliptical earth, not used.
453 igdstmpl(8)=kgds(2) ! octs 31-34; # "i" points
454 igdstmpl(9)=kgds(3) ! octs 35-38; # "j" points
455 igdstmpl(10)=1 ! octs 39-42; basic angle
456 igdstmpl(11)=10**6 ! octs 43-46; subdivisions of basic angle
457
458 scale=float(igdstmpl(10)*igdstmpl(11))
459
460 igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point
461
462 if (lon11 < 0) then
463 igdstmpl(13)=nint((lon11+360.)*scale) ! octs 51-54; lon of first grid point
464 else
465 igdstmpl(13)=nint(lon11*scale)
466 endif
467
468 igdstmpl(14) = 0 ! oct 55; resolution and component flags
469 if (btest(kgds(6),7)) igdstmpl(14) = 48
470 if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8
471
472 igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point
473
474 if (lonlast < 0) then
475 igdstmpl(16)=nint((lonlast+360.)*scale) ! octs 60-63; lon of last grid point
476 else
477 igdstmpl(16)=nint(lonlast*scale)
478 endif
479
480 igdstmpl(17)= nint(360.0/float(kgds(2)-1)*scale) ! octs 64-67; di of grid
481 igdstmpl(18)= kgds(3)/2 ! octs 68-71; # grid pts between pole and equator
482
483 igdstmpl(19)=0 ! oct 72; scanning mode flag
484 if(btest(kgds(11),7)) igdstmpl(19)=128
485 if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64
486 if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32
487
488 igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1
489 igds(2) = kgds(2)*kgds(3) ! num grid points
490 igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid)
491 igds(4) = 0 ! regular grid, no appended list
492 igds(5) = 40 ! gaussian
493
494 ngrdpts = igds(2)
495
496! These variables used for non-regular grids. We are using regular grids
497! (igds(3) equals 0).
498
499 idefnum=1
500 ideflist=0
501
502 elseif (kgds(1) == 203 .or. kgds(1) == 205) then
503
504 igdstmpl(1)=5 ! oct 15; shape of the earth, wgs84
505 igdstmpl(2)=255 ! oct 16; scale factor of radius of spherical earth, not used.
506 igdstmpl(3)=-1 ! octs 17-20; scale value of radius of spherical earth, not used.
507 igdstmpl(4)=255 ! oct 21; scale factor of major axis of elliptical earth, not used.
508 igdstmpl(5)=-1 ! octs 22-25; scaled value of major axis of elliptical earth, not used.
509 igdstmpl(6)=255 ! oct 26; scale factor of minor axis of elliptical earth, not used.
510 igdstmpl(7)=-1 ! octs 27-30; scaled value of minor axis of elliptical earth, not used.
511 igdstmpl(8)=kgds(2) ! octs 31-34; # "i" points
512 igdstmpl(9)=kgds(3) ! octs 35-38; # "j" points
513 igdstmpl(10)=1 ! octs 39-42; basic angle
514 igdstmpl(11)=10**6 ! octs 43-46; subdivisions of basic angle
515
516 scale=float(igdstmpl(10)*igdstmpl(11))
517
518 igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point
519
520 if (lon11 < 0) then
521 igdstmpl(13)=nint((lon11+360.)*scale) ! octs 51-54; lon of first grid point
522 else
523 igdstmpl(13)=nint(lon11*scale)
524 endif
525
526 igdstmpl(14) = 0 ! oct 55; resolution and component flags
527 if (btest(kgds(6),7)) igdstmpl(14) = 48
528 if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8
529
530 igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point
531
532 if (lonlast < 0) then
533 igdstmpl(16)=nint((lonlast+360.)*scale) ! octs 60-63; lon of last grid point
534 else
535 igdstmpl(16)=nint(lonlast*scale)
536 endif
537
538 if (kgds(1) == 203) igdstmpl(17)= nint(float(kgds(9))*scale/500.) ! octs 64-67; di of grid.
539 ! iplib "e" grid convention
540 ! is 1/2 the grib convention.
541 if (kgds(1) == 205) igdstmpl(17)= nint(float(kgds(9))*scale/1000.) ! octs 64-67; di of grid
542
543 igdstmpl(18)= nint(float(kgds(10))*scale/1000.) ! octs 68-71; dj of grid
544
545 igdstmpl(19)=0 ! oct 72; scanning mode flag
546 if(btest(kgds(11),7)) igdstmpl(19)=128
547 if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64
548 if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32
549 if (kgds(1) == 203) igdstmpl(19)=igdstmpl(19) + 4
550
551 igdstmpl(20) = nint(float(kgds(7)-90000)*scale/1000.) ! octs 73-76; lat of south pole of projection
552
553 if (kgds(8) < 0) then
554 igdstmpl(21) = nint(float(kgds(8)+360000)*scale/1000.) ! octs 77-80; long of southern pole of projection.
555 else
556 igdstmpl(21) = nint(float(kgds(8))*scale/1000.) ! octs 77-80; long of southern pole of projection.
557 endif
558
559 igdstmpl(22)=0 ! octs 81-84; angle of rotation of projection
560
561 igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1
562 igds(2) = kgds(2)*kgds(3) ! num grid points
563 igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid)
564 igds(4) = 0 ! regular grid, no appended list
565 igds(5) = 1 ! rotated lat/lon
566
567 ngrdpts = igds(2)
568
569! These variables used for non-regular grids. We are using regular grids
570! (igds(3) equals 0).
571
572 idefnum=1
573 ideflist=0
574
575 end if
576
577! Section 4 - product definition section
578
579 ipdsnum = 0 ! pds template number - table 4.0
580
581 ipdstmpl(1)= 1 ! oct 10; parameter category
582! note!! to use a parmeter number >= 192 you must set the local table to '1'
583 ipdstmpl(2)= 42 ! oct 11; parameter
584 ipdstmpl(3)= 0 ! oct 12; type of generating process
585 ipdstmpl(4)= 255 ! oct 13; background generating process identifier
586 ipdstmpl(5)= 84 ! oct 14; analysis generating process identifier
587 ipdstmpl(6)= 0 ! octs 15-16; hours after ob cutoff
588 ipdstmpl(7)= 0 ! oct 17; minutes after ob cutoff
589 ipdstmpl(8)= 1 ! oct 18; unit of time range
590 ipdstmpl(9)= 0 ! octs 19-22; forecast time in units defined by oct 18
591 ipdstmpl(10)=1 ! oct 23; type of first fixed surface
592 ipdstmpl(11)=0 ! oct 24; scale factor of first fixed surface
593 ipdstmpl(12)=0 ! octs 25-28; scale value of first fixed surface
594 ipdstmpl(13)=255 ! oct 29; type of second fixed surface
595 ipdstmpl(14)=255 ! oct 30; scale factor of second fixed surface
596 ipdstmpl(15)=-2147483647 ! octs 31-34; scaled value of second fixed surface
597 ! note! for these particular octets, using -1 as
598 ! missing does not work because -1 may be an actual
599 ! scaled value. after looking thru the g2 library
600 ! and some trial and error, i determined that missing
601 ! is minus 2**31-1.
602
603 end subroutine init_grib2
604
610 subroutine grib2_null(gfld)
611
612 use grib_mod
613
614 implicit none
615
616 type(gribfield), intent(inout) :: gfld
617
618 nullify(gfld%idsect)
619 nullify(gfld%local)
620 nullify(gfld%list_opt)
621 nullify(gfld%igdtmpl)
622 nullify(gfld%ipdtmpl)
623 nullify(gfld%coord_list)
624 nullify(gfld%idrtmpl)
625 nullify(gfld%bmap)
626 nullify(gfld%fld)
627
628 end subroutine grib2_null
629
635 subroutine grib2_free(gfld)
636
637 use grib_mod
638
639 implicit none
640
641 type(gribfield), intent(inout) :: gfld
642
643 if (associated(gfld%idsect)) deallocate(gfld%idsect)
644 if (associated(gfld%local)) deallocate(gfld%local)
645 if (associated(gfld%list_opt)) deallocate(gfld%list_opt)
646 if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl)
647 if (associated(gfld%ipdtmpl)) deallocate(gfld%ipdtmpl)
648 if (associated(gfld%coord_list)) deallocate(gfld%coord_list)
649 if (associated(gfld%idrtmpl)) deallocate(gfld%idrtmpl)
650 if (associated(gfld%bmap)) deallocate(gfld%bmap)
651 if (associated(gfld%fld)) deallocate(gfld%fld)
652
653 end subroutine grib2_free
grib2_check
subroutine grib2_check(kgds, igdstmplen)
Determine length of grib2 gds template array, which is a function of the map projection.
Definition grib_utils.F90:351
skgb2
subroutine skgb2(lugb, iseek, mseek, lskip, lgrib, i1)
Determine whether file is grib or not.
Definition grib_utils.F90:76
grib2_free
subroutine grib2_free(gfld)
Deallocate the grib2 gribfield pointers.
Definition grib_utils.F90:636
init_grib2
subroutine init_grib2(century, year, month, day, hour, kgds, lat11, latlast, lon11, lonlast, listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl, igdstmplen, idefnum, ideflist, ngrdpts)
Initialize grib2 arrays required by the ncep g2 library according to grib1 gds information.
Definition grib_utils.F90:401
gdt_to_gds
subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
Convert from the grib2 grid description template array used by the ncep grib2 library,...
Definition grib_utils.F90:140
grib2_null
subroutine grib2_null(gfld)
Nullify the grib2 gribfield pointers.
Definition grib_utils.F90:611
grib_check
subroutine grib_check(file_name, isgrib)
Determine whether file is grib or not.
Definition grib_utils.F90:24