| blob | 77da8e878430c7e0be437b670719bd51ded757fb |
1 !*****************************************************************************!
2 ! Subroutine RD_GRIB2 !
3 ! !
4 ! Purpose: !
5 ! Read one record from the input GRIB2 file. Based on the information in !
6 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
7 ! the GRIB2 record is one to process or to skip. If the field is one we !
8 ! want to keep, extract the data from the GRIB2 record, and store the data !
9 ! in the ungrib memory structure. !
10 ! !
11 ! Argument list: !
12 ! Input: !
13 ! junit : "Unit Number" to open and read from. Not really a Fortran !
14 ! unit number, since we do not do Fortran I/O for the GRIB2 !
15 ! files. Nor is it a UNIX File Descriptor returned from a C !
16 ! OPEN statement. It is really just an array index to the !
17 ! array (IUARR) where the UNIX File Descriptor values are !
18 ! stored. !
19 ! gribflnm : File name to open, if it is not already open. !
20 ! debug_level : Integer for various amounts of printout. !
21 ! !
22 ! Output: !
23 ! !
24 ! hdate : The (up to)19-character date of the field to process. !
25 ! grib_edition : Version of the gribfile (1 or 2) !
26 ! ireaderr : Error flag: 0 - no error on read from GRIB2 file. !
27 ! 1 - Hit the end of the GRIB2 file. !
28 ! 2 - The file GRIBFLNM we tried to open does !
29 ! not exist. !
30 ! !
31 ! !
32 ! Author: Paula McCaslin, NOAA/FSL, Sept 2004 !
33 ! Code is based on code developed by Steve Gilbert NCEP & Kevin Manning NCAR !
34 ! Adapted for WPS: Jim Bresch, NCAR/MMM. Sept 2006 !
35 !*****************************************************************************!
37 SUBROUTINE rd_grib2(junit, gribflnm, hdate,
38 & grib_edition, ireaderr, debug_level)
40 use grib_mod
41 use params
42 use table ! Included to define g2code
43 use gridinfo ! Included to define map%
44 use storage_module ! Included sub put_storage
45 use module_debug
47 real, allocatable, dimension(:) :: hold_array
48 parameter(msk1=32000,msk2=4000)
49 character(len=1),allocatable,dimension(:) :: cgrib
50 integer :: listsec0(3)
51 integer :: listsec1(13)
52 integer year, month, day, hour, minute, second, fcst
53 character(len=*) :: gribflnm
54 character(len=*) :: hdate
55 character(len=8) :: pabbrev
56 character(len=20) :: labbrev
57 character(len=80) :: tabbrev
58 integer :: lskip, lgrib
59 integer :: junit, itot, icount, iseek
60 integer :: grib_edition
61 integer :: i, j, ireaderr, ith , debug_level
62 integer :: currlen
63 logical :: unpack, expand
64 type(gribfield) :: gfld
65 ! For subroutine put_storage
66 real :: level
67 real :: scale_factor
68 integer :: iplvl
69 character (len=9) :: my_field
70 character (len=8) :: tmp8
71 ! For subroutine outout
72 integer , parameter :: maxlvl = 100
73 real , dimension(maxlvl) :: plvl
74 integer :: nlvl
75 integer , dimension(maxlvl) :: level_array
76 logical :: first = .true.
78 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 C SET ARGUMENTS
81 unpack=.true.
82 expand=.true.
83 hdate = '0000-00-00_00:00:00'
84 ierr=0
85 itot=0
86 icount=0
87 iseek=0
88 lskip=0
89 lgrib=0
90 currlen=0
91 ith=1
92 scale_factor = 1e6
93 call mprintf(.true.,DEBUG,"Begin rd_grib2", newline=.true.)
95 !/* IOS Return Codes from BACIO: */
96 !/* 0 All was well */
97 !/* -1 Tried to open read only _and_ write only */
98 !/* -2 Tried to read and write in the same call */
99 !/* -3 Internal failure in name processing */
100 !/* -4 Failure in opening file */
101 !/* -5 Tried to read on a write-only file */
102 !/* -6 Failed in read to find the 'start' location */
103 !/* -7 Tried to write to a read only file */
104 !/* -8 Failed in write to find the 'start' location */
105 !/* -9 Error in close */
106 !/* -10 Read or wrote fewer data than requested */
108 !if ireaderr =1 we have hit the end of a file.
109 !if ireaderr =2 we have hit the end of all the files.
112 ! Open a byte-addressable file.
113 CALL BAOPENR(junit,gribflnm,IOS)
114 first = .true.
115 if (ios.eq.0) then
116 VERSION: do
118 ! Search opend file for the next GRIB2 messege (record).
119 call skgb(junit,iseek,msk1,lskip,lgrib)
121 ! Check for EOF, or problem
122 if (lgrib.eq.0) then
123 exit
124 endif
126 ! Check size, if needed allocate more memory.
127 if (lgrib.gt.currlen) then
128 if (allocated(cgrib)) deallocate(cgrib)
129 allocate(cgrib(lgrib),stat=is)
130 !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
131 currlen=lgrib
132 endif
134 ! Read a given number of bytes from unblocked file.
135 call baread(junit,lskip,lgrib,lengrib,cgrib)
137 call mprintf ((lgrib.ne.lengrib),ERROR,
138 & "rd_grib2: IO Error. %i .ne. %i ", newline=.true.,
139 & i1=lgrib,i2=lengrib)
141 iseek=lskip+lgrib
142 icount=icount+1
144 call mprintf (.true.,DEBUG,
145 & "G2 GRIB MESSAGE %i starts at %i ", newline=.true.,
146 & i1=icount,i2=lskip+1)
148 ! Unpack GRIB2 field
149 call gb_info(cgrib,lengrib,listsec0,listsec1,
150 & numfields,numlocal,maxlocal,ierr)
151 call mprintf((ierr.ne.0),ERROR,
152 & " ERROR querying GRIB2 message = %i",newline=.true.,i1=ierr)
153 itot=itot+numfields
155 grib_edition=listsec0(2)
156 if (grib_edition.ne.2) then
157 exit VERSION
158 endif
160 ! Additional print statments for developer.
161 !MGD if ( debug_level .GT. 100 ) then
162 !MGD print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
163 !MGD print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
164 !MGD print *,'G2 Contains ',numlocal,' Local Sections ',
165 !MGD & ' and ',numfields,' data fields.'
166 !MGD endif
169 ! ----
170 ! Once per file fill in date, model and projection values.
172 if (first) then
173 first = .false.
175 ! Build the 19-character date string, based on GRIB2 header date
176 ! and time information, including forecast time information:
178 n=1
179 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
180 year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
181 month =gfld%idsect(7) ! MONTH OF THE DATA
182 day =gfld%idsect(8) ! DAY OF THE DATA
183 hour =gfld%idsect(9) ! HOUR OF THE DATA
184 minute=gfld%idsect(10) ! MINUTE OF THE DATA
185 second=gfld%idsect(11) ! SECOND OF THE DATA
187 fcst = 0
189 ! Extract forecast time.
190 if ( gfld%ipdtmpl(8) .eq. 1 ) then ! time units are hours
191 fcst = gfld%ipdtmpl(9)
192 else if ( gfld%ipdtmpl(8) .eq. 0 ) then ! minutes
193 fcst = gfld%ipdtmpl(9) / 60.
194 else if ( gfld%ipdtmpl(8) .eq. 2 ) then ! days
195 fcst = gfld%ipdtmpl(9) * 24.
196 else
197 call mprintf(.true.,ERROR,
198 & "Time unit in ipdtmpl(8), %i is not suported",
199 & newline=.true.,i1=gfld%ipdtmpl(8))
200 endif
202 ! Compute valid time.
204 !print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
205 !print *, 'hhmm ',gfld%idsect(9),gfld%idsect(10)
207 call build_hdate(hdate,year,month,day,hour,minute,second)
208 call mprintf(.true.,DEBUG,"G2 hdate = %s ", newline=.true.,
209 & s1=hdate)
210 call geth_newdate(hdate,hdate,3600*fcst)
211 call mprintf(.true.,DEBUG,"G2 hdate (fcst?) = %s ",
212 & newline=.true., s1=hdate)
214 !--
216 ! Indicator of the source (center) of the data.
217 icenter = gfld%idsect(1)
219 ! Indicator of model (or whatever) which generated the data.
220 iprocess = gfld%ipdtmpl(5)
223 if (icenter.eq.7) then
224 if (iprocess.eq.83 .or. iprocess.eq.84) then
225 map%source = 'NCEP MESO NAM Model'
226 elseif (iprocess.eq.81) then
227 map%source = 'NCEP GFS Analysis'
228 elseif (iprocess.eq.82) then
229 map%source = 'NCEP GFS GDAS/FNL'
230 elseif (iprocess.eq.89) then
231 map%source = 'NCEP NMM '
232 elseif (iprocess.eq.96) then
233 map%source = 'NCEP GFS Model'
234 elseif (iprocess.eq.86 .or. iprocess.eq.100) then
235 map%source = 'NCEP RUC Model' ! 60 km
236 elseif (iprocess.eq.101) then
237 map%source = 'NCEP RUC Model' ! 40 km
238 elseif (iprocess.eq.105) then
239 map%source = 'NCEP RUC Model' ! 20 km
240 elseif (iprocess.eq.107) then
241 map%source = 'NCEP GEFS'
242 elseif (iprocess.eq.109) then
243 map%source = 'NCEP RTMA'
244 elseif (iprocess.eq.140) then
245 map%source = 'NCEP NARR'
246 elseif (iprocess.eq.44) then
247 map%source = 'NCEP SST Analysis'
248 elseif (iprocess.eq.70) then
249 map%source = 'GFDL Hurricane Model'
250 elseif (iprocess.eq.80) then
251 map%source = 'NCEP GFS Ensemble'
252 elseif (iprocess.eq.107) then ! renumbered as of 23 Feb 2010
253 map%source = 'NCEP GFS Ensemble'
254 elseif (iprocess.eq.129) then
255 map%source = 'NCEP GODAS'
256 elseif (iprocess.eq.25) then
257 map%source = 'NCEP SNOW COVER ANALYSIS'
258 else
259 map%source = 'unknown model from NCEP'
260 call mprintf(.true.,STDOUT,
261 & "unknown model from NCEP %i ",newline=.true.,
262 & i1=iprocess)
263 call mprintf(.true.,LOGFILE,
264 & "unknown model from NCEP %i ",newline=.true.,
265 & i1=iprocess)
266 end if
267 else if (icenter .eq. 57) then
268 if (iprocess .eq. 87) then
269 map%source = 'AFWA AGRMET'
270 else
271 map%source = 'AFWA'
272 endif
273 else if ( icenter .eq. 58 ) then
274 map%source = 'US Navy FNOC'
275 else if (icenter .eq. 59) then
276 if (iprocess.eq.120) then
277 map%source = 'NOAA GSD Rapid Refresh'
278 else
279 map%source = 'NOAA GSD'
280 endif
281 else if (icenter .eq. 60) then
282 map%source = 'NCAR'
283 else if (icenter .eq. 98) then
284 map%source = 'ECMWF'
285 else if (icenter .eq. 34) then
286 map%source = 'JMA'
287 else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
288 map%source = 'UKMO'
289 else
290 map%source = 'unknown model and orig center'
291 end if
292 call mprintf(.true.,DEBUG,"G2 source is = %s ",
293 & newline=.true., s1=map%source)
295 !--
297 ! Store information about the grid containing the data.
298 ! This stuff gets stored in the MAP variable, as defined in
299 ! module GRIDINFO.
301 map%startloc = 'SWCORNER'
302 map%grid_wind = .true.
304 if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
305 map%igrid = 0
306 map%nx = gfld%igdtmpl(8)
307 map%ny = gfld%igdtmpl(9)
308 map%dx = gfld%igdtmpl(17)
309 map%dy = gfld%igdtmpl(18)
310 map%lat1 = gfld%igdtmpl(12)
311 map%lon1 = gfld%igdtmpl(13)
312 write(tmp8,'(b8.8)') gfld%igdtmpl(14)
313 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
314 map%r_earth = earth_radius (gfld%igdtmpl(1),
315 & gfld%igdtmpl(2),gfld%igdtmpl(3))
317 ! Fix for NCEP 1/12 degree grids (e.g. rtgsst)
318 if (icenter .eq. 7 .and. map%dx .eq. 83000. .and. map%nx
319 & .eq. 4320) then
320 map%lat1 = 89958333.
321 map%lon1 = 41667.
322 map%dx = 83333.333 * sign(1.0,map%dx)
323 map%dy = 83333.333 * sign(1.0,map%dy)
324 endif
326 if ((gfld%igdtmpl(10) .eq. 0).OR.
327 & (gfld%igdtmpl(10) .eq. 255)) THEN
328 ! Scale lat/lon values to 0-180, default range is 1e6.
329 map%lat1 = map%lat1/scale_factor
330 map%lon1 = map%lon1/scale_factor
331 ! Scale dx/dy values to degrees, default range is 1e6.
332 map%dx = map%dx/scale_factor
333 map%dy = map%dy/scale_factor
334 else
335 ! Basic angle and subdivisions are non-zero (not tested)
336 map%lat1 = map%lat1 *
337 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
338 map%lon1 = map%lon1 *
339 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
340 map%dx = map%dx *
341 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
342 map%dy = map%dy *
343 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
344 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
345 &has not been tested, continuing anyway")
346 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
347 & has not been tested, continuing anyway")
348 endif
351 ! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
352 ! In WPS, the sign of dy indicates the direction of the scan.
353 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
354 read(tmp8,'(1x,i1)') jscan
355 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
356 map%dy = -1. * map%dy
357 endif
358 ! if ( map%lat1 .gt. gfld%igdtmpl(15) .and.
359 ! & map%dy .gt. 0. ) then
360 ! map%dy = -1. * map%dy
361 ! write(6,*) 'Resetting map%dy for iprocess = ',iprocess
362 ! endif
364 elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
365 map%igrid = 5
366 map%nx = gfld%igdtmpl(8)
367 map%ny = gfld%igdtmpl(9)
368 map%lov = gfld%igdtmpl(14) / scale_factor
369 map%truelat1 = 60.
370 map%truelat2 = 91.
371 map%dx = gfld%igdtmpl(15) / scale_factor
372 map%dy = gfld%igdtmpl(16) / scale_factor
373 map%lat1 = gfld%igdtmpl(10) / scale_factor
374 map%lon1 = gfld%igdtmpl(11) / scale_factor
375 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
376 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
377 map%r_earth = earth_radius (gfld%igdtmpl(1),
378 & gfld%igdtmpl(2),gfld%igdtmpl(3))
380 elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
381 map%igrid = 3
382 map%nx = gfld%igdtmpl(8)
383 map%ny = gfld%igdtmpl(9)
384 map%lov = gfld%igdtmpl(14) / scale_factor
385 map%truelat1 = gfld%igdtmpl(19) / scale_factor
386 map%truelat2 = gfld%igdtmpl(20) / scale_factor
387 map%dx = gfld%igdtmpl(15) / scale_factor
388 map%dy = gfld%igdtmpl(16) / scale_factor
389 map%lat1 = gfld%igdtmpl(10) / scale_factor
390 map%lon1 = gfld%igdtmpl(11) / scale_factor
391 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
392 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
393 map%r_earth = earth_radius (gfld%igdtmpl(1),
394 & gfld%igdtmpl(2),gfld%igdtmpl(3))
396 elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
397 map%igrid = 4
398 map%nx = gfld%igdtmpl(8) ! Ni - # of points along a parallel
399 map%ny = gfld%igdtmpl(9) ! Nj - # of points along meridian
400 map%dx = gfld%igdtmpl(17) ! Di - i direction increment
401 map%dy = gfld%igdtmpl(18) ! N - # of parallels between pole and equator
402 map%lat1 = gfld%igdtmpl(12) ! La1 - lat of 1st grid point
403 map%lon1 = gfld%igdtmpl(13) ! Lo1 - lon of 1st grid point
404 write(tmp8,'(b8.8)') gfld%igdtmpl(14) ! resolution/component flag
405 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
406 map%r_earth = earth_radius (gfld%igdtmpl(1),
407 & gfld%igdtmpl(2),gfld%igdtmpl(3))
409 ! Scale dx/dy values to degrees, default range is 1e6.
410 if (map%dx.gt.10000) then
411 map%dx = map%dx/scale_factor
412 endif
413 if (map%dy.gt.10000) then
414 map%dy = (map%dy/scale_factor)*(-1)
415 endif
417 ! Scale lat/lon values to 0-180, default range is 1e6.
418 if (map%lat1.ge.scale_factor) then
419 map%lat1 = map%lat1/scale_factor
420 endif
421 if (map%lon1.ge.scale_factor) then
422 map%lon1 = map%lon1/scale_factor
423 endif
424 if ( debug_level .gt. 2 ) then
425 call mprintf(.true.,DEBUG,
426 & "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
427 & newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
428 & i1=nint(map%dy))
429 end if
431 else
432 call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
433 & newline=.true.,i1=gfld%igdtnum)
434 call mprintf(.true.,STDOUT,
435 & "see Code Table 3.1: Grid Definition Template Number",
436 & newline=.true.)
437 call mprintf(.true.,LOGFILE,
438 & "GRIB2 Unknown Projection: %i",
439 & newline=.true.,i1=gfld%igdtnum)
440 call mprintf(.true.,LOGFILE,
441 & "see Code Table 3.1: Grid Definition Template Number",
442 & newline=.true.)
443 endif
445 call mprintf(.true.,DEBUG,"G2 igrid = %i , dx = %f , dy = %
446 &f ", newline=.true., i1 = map%igrid, f1=map%dx, f2=map%dy)
448 if (icenter.eq.7) then
449 call ncep_grid_num (gfld%igdtnum)
450 endif
452 ! Deallocate arrays decoding GRIB2 record.
453 call gf_free(gfld)
454 endif
456 ! ----
458 ! Continue to unpack GRIB2 field.
459 do n=1,numfields ! e.g. U and V would =2, otherwise its usually =1
460 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
461 if (ierr.ne.0) then
462 write(*,*) ' ERROR extracting field gf_getfld = ',ierr
463 cycle
464 endif
466 ! The JMA GSM has two different grids in the same GRIB file, so we need
467 ! to process the map info for each field separately. If any other centers do
468 ! this, then processing will need to be added here, too.
470 if (icenter .eq. 34 .and. gfld%igdtnum.eq.0) then
471 map%nx = gfld%igdtmpl(8)
472 map%ny = gfld%igdtmpl(9)
473 map%dx = gfld%igdtmpl(17)
474 map%dy = gfld%igdtmpl(18)
475 ! Scale dx/dy values to degrees, default range is 1e6.
476 if (map%dx.gt.10000) then
477 map%dx = map%dx/scale_factor
478 endif
479 if (map%dy.gt.10000) then
480 map%dy = map%dy/scale_factor
481 endif
482 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
483 read(tmp8,'(1x,i1)') jscan
484 write(0,*) 'gfld%igdtmpl(19) = ',gfld%igdtmpl(19),
485 & ' jscan = ',jscan
486 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
487 map%dy = -1. * map%dy
488 endif
489 endif ! JMA spectral
491 ! ------------------------------------
492 ! Additional print information for developer.
493 if ( debug_level .GT. 1000 ) then
494 !MGD print *
495 !MGD print *,'G2 FIELD ',n
496 !MGD if (n==1) then
497 !MGD print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
498 !MGD print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
499 !MGD endif
500 !MGD if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
501 !MGD print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
502 !MGD endif
503 !MGD print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
504 !MGD & gfld%numoct_opt,gfld%interp_opt,
505 !MGD & gfld%igdtnum
506 !MGD print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
507 !MGD & (gfld%igdtmpl(j),j=1,gfld%igdtlen)
508 !MGD if ( gfld%num_opt .eq. 0 ) then
509 !MGD print *,'G2 NO Section 3 List Defining No. of Data Points.'
510 !MGD else
511 !MGD print *,'G2 Section 3 Optional List: ',
512 !MGD & (gfld%list_opt(j),j=1,gfld%num_opt)
513 !MGD endif
514 !MGD print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
515 !MGD & (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
517 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
518 & gfld%ipdtmpl(2))
519 !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
520 !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
521 !MGD print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
523 !MGD if ( gfld%num_coord .eq. 0 ) then
524 !MGD print *,'G2 NO Optional Vertical Coordinate List.'
525 !MGD else
526 !MGD print *,'G2 Section 4 Optional Coordinates: ',
527 !MGD & (gfld%coord_list(j),j=1,gfld%num_coord)
528 !MGD endif
529 if ( gfld%ibmap .ne. 255 ) then
530 call mprintf(.true.,DEBUG,
531 & 'G2 Num. of Data Points = %i with BIT-MAP %i',
532 & newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
533 else
534 call mprintf(.true.,DEBUG,
535 & 'G2 Num. of Data Points = %i NO BIT-MAP',
536 & newline=.true., i1=gfld%ndpts)
537 endif
538 !MGD print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
539 !MGD & (gfld%idrtmpl(j),j=1,gfld%idrtlen)
540 fldmax=gfld%fld(1)
541 fldmin=gfld%fld(1)
542 sum=gfld%fld(1)
543 do j=2,gfld%ndpts
544 if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
545 if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
546 sum=sum+gfld%fld(j)
547 enddo ! gfld%ndpts
549 !MGD print *,'G2 Data Values:'
550 call mprintf(.true.,DEBUG,'G2 MIN=%f AVE=%f MAX=%f',
551 & newline=.true., f1=fldmin, f2=sum/gfld%ndpts, f3=fldmax)
552 !do j=1,gfld%ndpts\20
553 ! write(*,*) j, gfld%fld(j)
554 !enddo
555 endif ! Additional Print information
556 ! ------------------------------------
558 ! do i = 1, maxvar
559 ! write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
560 ! enddo
561 !MGD if (debug_level .gt. 50) then
562 !MGD write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
563 !MGD & gfld%ipdtmpl(2),gfld%ipdtmpl(10)
564 !MGD endif
565 call mprintf(.true.,DEBUG,"G2 Searching the g2code array (Vta
566 &ble) for this grib field %i %i %i %i ", newline=.true.,
567 & i1 = gfld%discipline, i2 = gfld%ipdtmpl(1),
568 & i3 = gfld%ipdtmpl(2), i4 = gfld%ipdtmpl(10) )
570 ! Test this data record against list of desired variables
571 ! found in Vtable.
572 ! ----
573 MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
574 ! maxvar is defined in table.mod
576 if ( gfld%discipline .eq. g2code(1,i) .and. !Discipline
577 & gfld%ipdtmpl(1) .eq. g2code(2,i) .and. !Category
578 & gfld%ipdtmpl(2) .eq. g2code(3,i) .and. !Parameter
579 & gfld%ipdtmpl(10) .eq. g2code(4,i) .and. !Elevation
580 & gfld%ipdtnum .eq. g2code(5,i)) then !Template
582 call gf_free(gfld)
583 call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
584 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
585 & gfld%ipdtmpl(2))
587 !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
588 my_field=namvar(i)
590 !MGD if (debug_level .gt. 50) then
591 !MGD write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
592 !MGD write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
593 !MGD endif
596 ! need to match up soil levels with those requested.
597 ! For the Vtable levels, -88 = all levels, -99 = missing. The units
598 ! vary depending on the level code (e.g. 106 = cm, 103 = m).
599 if ( gfld%ipdtmpl(10) .eq. 106 ) then
600 TMP8LOOP: do j = 1, maxvar
601 if ((g2code(4,j) .eq. 106) .and.
602 & (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
603 & (gfld%ipdtmpl(12) .eq. level1(j)) .and.
604 & ((gfld%ipdtmpl(15) .eq. level2(j)) .or.
605 & (level2(j) .le. -88))) then
606 my_field = namvar(j)
607 exit TMP8LOOP
608 endif
609 enddo TMP8LOOP
610 if (j .gt. maxvar ) then
611 write(6,'(a,i6,a,i6,a)') 'Subsoil level ',
612 & gfld%ipdtmpl(12),' to ',gfld%ipdtmpl(15),
613 & ' in the GRIB2 file, was not found in the Vtable'
614 cycle MATCH_LOOP
615 endif
616 !MGD if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
617 endif
619 ! Level (eg. 10000 mb)
620 if(gfld%ipdtmpl(10).eq.100) then
621 ! Pressure level (range from 1000mb to 0mb)
622 level=gfld%ipdtmpl(12) *
623 & (10. ** (-1. * gfld%ipdtmpl(11)))
624 elseif(gfld%ipdtmpl(10).eq.105) then
625 ! Hybrid level (range from 1 to N)
626 level=gfld%ipdtmpl(12)
627 elseif(gfld%ipdtmpl(10).eq.104) then
628 ! Sigma level (range from 10000 to 0)
629 level=gfld%ipdtmpl(12)
630 elseif(gfld%ipdtmpl(10).eq.101) then
631 ! MSL
632 level=201300.
633 elseif(gfld%ipdtmpl(10).eq.106.or.
634 & gfld%ipdtmpl(10).eq.1) then
635 ! Misc near ground/surface levels
636 level=200100.
637 elseif(gfld%ipdtmpl(10).eq.6) then
638 ! Level of Max wind
639 level=6. ! .06 mb should never be used by anyone
640 elseif(gfld%ipdtmpl(10).eq.7) then
641 ! Tropopause
642 level=7. ! .07 mb should never be used by anyone
643 else
644 ! Misc near ground/surface levels
645 level=200100.
646 endif
647 iplvl = int(level)
649 ! Store the unpacked 2D slab from the GRIB2 record
650 allocate(hold_array(gfld%ngrdpts))
651 do j=1,gfld%ngrdpts
652 hold_array(j)=gfld%fld(j)
653 enddo
655 ! Some grids need to be reordered. Until we get an example, this is
656 ! a placeholder
657 ! call reorder_it (hold_array, map%nx, map%ny, map%dx,
658 ! & map%dy, iorder)
660 ! When we have reached this point, we have a data array ARRAY
661 ! which has some data we want to save, with field name FIELD
662 ! at pressure level LEVEL (Pa). Dimensions of this data are
663 ! map%nx and map%ny. Put this data into storage.
665 !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
666 !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
667 call put_storage(iplvl,my_field,
668 & reshape(hold_array(1:map%nx*map%ny),
669 & (/map%nx, map%ny/)), map%nx,map%ny)
670 deallocate(hold_array)
672 ! If Specific Humidity is present on hybrid levels AND
673 ! upper-air RH is missing, see if we can compute RH from
674 ! Specific Humidity.
675 if (.not. is_there(iplvl, 'RH') .and.
676 & is_there(iplvl, 'SH') .and.
677 & is_there(iplvl, 'TT') .and.
678 & is_there(iplvl, 'P')) then
679 call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
680 !call llstor_remove(iplvl, 'SH') !We are done with SH
681 endif
683 ith=ith+1
684 exit MATCH_LOOP
686 endif ! Selected param.
689 enddo MATCH_LOOP
691 ! Deallocate arrays decoding GRIB2 record.
692 call gf_free(gfld)
694 enddo ! 1,numfields
697 enddo VERSION ! skgb
700 if ( debug_level .gt. 100 ) then
701 call mprintf (.true.,DEBUG,
702 & "G2 total number of fields found = %i ",newline=.true.,i1=itot)
703 end if
705 CALL BACLOSE(junit,IOS)
707 nullify(gfld%local) ! must be nullified before opening next file
708 ireaderr=1
709 else
710 call mprintf (.true.,DEBUG,"open status failed because %i ",
711 & newline=.true., i1=ios)
712 hdate = '9999-99-99_99:99:99'
713 ireaderr=2
714 endif ! ireaderr check
716 END subroutine rd_grib2
718 !*****************************************************************************!
719 ! Subroutine edition_num !
720 ! !
721 ! Purpose: !
722 ! Read one record from the input GRIB2 file. Based on the information in !
723 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
724 ! the GRIB2 record is one to process or to skip. If the field is one we !
725 ! want to keep, extract the data from the GRIB2 record, and pass the data !
726 ! back to the calling routine. !
727 ! !
728 ! Argument list: !
729 ! Input: !
730 ! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
731 ! unit number, since we do not do Fortran I/O for the GRIB2 !
732 ! files. Nor is it a UNIX File Descriptor returned from a C !
733 ! OPEN statement. It is really just an array index to the !
734 ! array (IUARR) where the UNIX File Descriptor values are !
735 ! stored. !
736 ! GRIB2FILE: File name to open, if it is not already open. !
737 ! !
738 ! Output: !
739 ! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
740 ! IERR : Error flag: 0 - no error on read from GRIB2 file. !
741 ! 1 - Hit the end of the GRIB2 file. !
742 ! 2 - The file GRIBFLNM we tried to open does !
743 ! not exist. !
744 ! Author: Paula McCaslin !
745 ! NOAA/FSL !
746 ! Sept 2004 !
747 !*****************************************************************************!
749 SUBROUTINE edition_num(junit, gribflnm,
750 & grib_edition, ireaderr)
752 use grib_mod
753 use params
754 use module_debug
756 parameter(msk1=32000,msk2=4000)
757 character(len=1),allocatable,dimension(:) :: cgrib
758 integer :: listsec0(3)
759 integer :: listsec1(13)
760 character(len=*) :: gribflnm
761 integer :: lskip, lgrib
762 integer :: junit
763 integer :: grib_edition
764 integer :: i, j, ireaderr
765 integer :: currlen
767 character(len=4) :: ctemp
768 character(len=4),parameter :: grib='GRIB',c7777='7777'
770 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
771 C SET ARGUMENTS
773 itot=0
774 icount=0
775 iseek=0
776 lskip=0
777 lgrib=0
778 currlen=0
780 !/* IOS Return Codes from BACIO: */
781 !/* 0 All was well */
782 !/* -1 Tried to open read only _and_ write only */
783 !/* -2 Tried to read and write in the same call */
784 !/* -3 Internal failure in name processing */
785 !/* -4 Failure in opening file */
786 !/* -5 Tried to read on a write-only file */
787 !/* -6 Failed in read to find the 'start' location */
788 !/* -7 Tried to write to a read only file */
789 !/* -8 Failed in write to find the 'start' location */
790 !/* -9 Error in close */
791 !/* -10 Read or wrote fewer data than requested */
793 !if ireaderr =1 we have hit the end of a file.
794 !if ireaderr =2 we have hit the end of all the files.
795 !if ireaderr =3 beginning characters 'GRIB' not found
797 ! Open a byte-addressable file.
798 CALL BAOPENR(junit,gribflnm,IOS)
799 if (ios.eq.0) then
801 ! Search opend file for the next GRIB2 messege (record).
802 call skgb(junit,iseek,msk1,lskip,lgrib)
804 ! Check for EOF, or problem
805 call mprintf((lgrib.eq.0),ERROR,
806 & "Grib2 file or date problem, stopping in edition_num.",
807 & newline=.true.)
809 ! Check size, if needed allocate more memory.
810 if (lgrib.gt.currlen) then
811 if (allocated(cgrib)) deallocate(cgrib)
812 allocate(cgrib(lgrib),stat=is)
813 currlen=lgrib
814 endif
816 ! Read a given number of bytes from unblocked file.
817 call baread(junit,lskip,lgrib,lengrib,cgrib)
819 ! Check for beginning of GRIB message in the first 100 bytes
820 istart=0
821 do j=1,100
822 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
823 if (ctemp.eq.grib ) then
824 istart=j
825 exit
826 endif
827 enddo
828 if (istart.eq.0) then
829 ireaderr=3
830 print*, "The beginning 4 characters >GRIB< were not found."
831 endif
833 ! Unpack Section 0 - Indicator Section to extract GRIB edition field
834 iofst=8*(istart+5)
835 call gbyte(cgrib,discipline,iofst,8) ! Discipline
836 iofst=iofst+8
837 call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
839 print *, 'ungrib - grib edition num', grib_edition
840 CALL BACLOSE(junit,IOS)
841 ireaderr=1
842 else if (ios .eq. -4) then
843 call mprintf(.true.,ERROR,
844 & "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
845 else
846 print *,'edition_num: open status failed because',ios,gribflnm
847 ireaderr=2
848 endif ! ireaderr check
850 END subroutine edition_num
852 !*****************************************************************************!
854 SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
855 ! Compute relative humidity from specific humidity in the upper air.
856 use storage_module
857 implicit none
858 integer :: ix, jx
859 integer :: iiplvl
860 real :: lat1, lon1, dx, dy
861 real, dimension(ix,jx) :: T, P, RH, Q
863 real, parameter :: svp1=611.2
864 real, parameter :: svp2=17.67
865 real, parameter :: svp3=29.65
866 real, parameter :: svpt0=273.15
867 real, parameter :: eps = 0.622
869 real startlat, startlon, deltalat, deltalon
871 call get_storage(iiplvl, 'P', P, ix, jx)
872 call get_storage(iiplvl, 'TT', T, ix, jx)
873 call get_storage(iiplvl, 'SH', Q, ix, jx)
875 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
877 call put_storage(iiplvl, 'RH', rh, ix, jx)
879 end subroutine g2_compute_rh_spechumd_upa
881 !*****************************************************************************!
883 subroutine ncep_grid_num (pnum)
884 !
885 ! Grib2 doesn't have a grid-number entry, so we have to figure it out
886 !
887 use gridinfo ! Included to define map%
888 integer :: pnum
889 real, parameter :: eps = .01
890 character (len=8) :: tmp8
892 ! write(6,*) 'begin ncep_grid_num'
893 ! write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
894 tmp8 = ' '
895 if (pnum .eq. 30) then
896 if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
897 write(tmp8,'("GRID 218")')
898 else if (abs(map%dx - 40.63525) .lt. eps
899 & .and. map%nx .eq. 185) then
900 write(tmp8,'("GRID 212")')
901 else if (abs(map%dx - 40.63525) .lt. eps
902 & .and. map%nx .eq. 151) then
903 write(tmp8,'("GRID 236")')
904 else if (abs(map%dx - 81.2705) .lt. eps
905 & .and. map%nx .eq. 93) then
906 write(tmp8,'("GRID 211")')
907 else if (abs (map%dx - 32.46341) .lt. eps
908 & .and. map%nx .eq. 349) then
909 write(tmp8,'("GRID 221")')
910 else if (abs(map%dx - 20.317625) .lt. eps
911 & .and. map%nx .eq. 301) then
912 write(tmp8,'("GRID 252")')
913 endif
914 else if (pnum .eq. 20) then
915 if (abs(map%dx - 15.0) .lt. eps) then
916 write(tmp8,'("GRID 88")')
917 endif
918 else if (pnum .eq. 0) then
919 if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
920 write(tmp8,'("GRID 3")')
921 else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
922 write(tmp8,'("GRID 4")')
923 endif
924 endif
925 map%source(25:32) = tmp8
926 ! write(6,*) 'map%source = ',map%source
927 end subroutine ncep_grid_num
928 !*****************************************************************************!
930 function earth_radius (icode, iscale, irad_m)
931 ! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
932 use module_debug
933 real :: earth_radius
934 integer :: icode
935 integer :: iscale, irad_m
936 if ( icode .eq. 0 ) then
937 earth_radius = 6367470. * .001
938 else if ( icode .eq. 1) then
939 earth_radius = 0.001 * float(irad_m) / 10**iscale
940 else if ( icode .eq. 6 ) then
941 earth_radius = 6371229. * .001
942 else if ( icode .eq. 8 ) then
943 earth_radius = 6371200. * .001
944 else
945 call mprintf(.true.,ERROR,
946 & "unknown earth radius for code %i",newline=.true.,i1=icode)
947 endif
948 end function earth_radius