| blob | f8a7956c68cfd76971e2abdd01f0118898377739 |
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. 125) then
277 map%source = 'NOAA GSD Rapid Refresh'
278 else if (iprocess .eq. 105) then
279 map%source = 'NOAA GSD'
280 else
281 print *,'Unknown GSD source'
282 stop
283 endif
284 else if (icenter .eq. 60) then
285 map%source = 'NCAR'
286 else if (icenter .eq. 98) then
287 map%source = 'ECMWF'
288 else if (icenter .eq. 34) then
289 map%source = 'JMA'
290 else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
291 map%source = 'UKMO'
292 else
293 map%source = 'unknown model and orig center'
294 end if
295 call mprintf(.true.,DEBUG,"G2 source is = %s ",
296 & newline=.true., s1=map%source)
298 !--
300 ! Store information about the grid containing the data.
301 ! This stuff gets stored in the MAP variable, as defined in
302 ! module GRIDINFO.
304 map%startloc = 'SWCORNER'
305 map%grid_wind = .true.
307 if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
308 map%igrid = 0
309 map%nx = gfld%igdtmpl(8)
310 map%ny = gfld%igdtmpl(9)
311 map%dx = gfld%igdtmpl(17)
312 map%dy = gfld%igdtmpl(18)
313 map%lat1 = gfld%igdtmpl(12)
314 map%lon1 = gfld%igdtmpl(13)
315 write(tmp8,'(b8.8)') gfld%igdtmpl(14)
316 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
317 map%r_earth = earth_radius (gfld%igdtmpl(1),
318 & gfld%igdtmpl(2),gfld%igdtmpl(3))
320 ! Fix for NCEP 1/12 degree grids (e.g. rtgsst)
321 if (icenter .eq. 7 .and. map%dx .eq. 83000. .and. map%nx
322 & .eq. 4320) then
323 map%lat1 = 89958333.
324 map%lon1 = 41667.
325 map%dx = 83333.333 * sign(1.0,map%dx)
326 map%dy = 83333.333 * sign(1.0,map%dy)
327 endif
329 if ((gfld%igdtmpl(10) .eq. 0).OR.
330 & (gfld%igdtmpl(10) .eq. 255)) THEN
331 ! Scale lat/lon values to 0-180, default range is 1e6.
332 map%lat1 = map%lat1/scale_factor
333 map%lon1 = map%lon1/scale_factor
334 ! Scale dx/dy values to degrees, default range is 1e6.
335 map%dx = map%dx/scale_factor
336 map%dy = map%dy/scale_factor
337 else
338 ! Basic angle and subdivisions are non-zero (not tested)
339 map%lat1 = map%lat1 *
340 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
341 map%lon1 = map%lon1 *
342 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
343 map%dx = map%dx *
344 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
345 map%dy = map%dy *
346 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
347 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
348 &has not been tested, continuing anyway")
349 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
350 & has not been tested, continuing anyway")
351 endif
354 ! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
355 ! In WPS, the sign of dy indicates the direction of the scan.
356 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
357 read(tmp8,'(1x,i1)') jscan
358 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
359 map%dy = -1. * map%dy
360 endif
361 ! if ( map%lat1 .gt. gfld%igdtmpl(15) .and.
362 ! & map%dy .gt. 0. ) then
363 ! map%dy = -1. * map%dy
364 ! write(6,*) 'Resetting map%dy for iprocess = ',iprocess
365 ! endif
367 elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
368 map%igrid = 5
369 map%nx = gfld%igdtmpl(8)
370 map%ny = gfld%igdtmpl(9)
371 map%lov = gfld%igdtmpl(14) / scale_factor
372 map%truelat1 = 60.
373 map%truelat2 = 91.
374 map%dx = gfld%igdtmpl(15) / scale_factor
375 map%dy = gfld%igdtmpl(16) / scale_factor
376 map%lat1 = gfld%igdtmpl(10) / scale_factor
377 map%lon1 = gfld%igdtmpl(11) / scale_factor
378 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
379 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
380 map%r_earth = earth_radius (gfld%igdtmpl(1),
381 & gfld%igdtmpl(2),gfld%igdtmpl(3))
383 elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
384 map%igrid = 3
385 map%nx = gfld%igdtmpl(8)
386 map%ny = gfld%igdtmpl(9)
387 map%lov = gfld%igdtmpl(14) / scale_factor
388 map%truelat1 = gfld%igdtmpl(19) / scale_factor
389 map%truelat2 = gfld%igdtmpl(20) / scale_factor
390 map%dx = gfld%igdtmpl(15) / scale_factor
391 map%dy = gfld%igdtmpl(16) / scale_factor
392 map%lat1 = gfld%igdtmpl(10) / scale_factor
393 map%lon1 = gfld%igdtmpl(11) / scale_factor
394 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
395 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
396 map%r_earth = earth_radius (gfld%igdtmpl(1),
397 & gfld%igdtmpl(2),gfld%igdtmpl(3))
399 elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
400 map%igrid = 4
401 map%nx = gfld%igdtmpl(8) ! Ni - # of points along a parallel
402 map%ny = gfld%igdtmpl(9) ! Nj - # of points along meridian
403 map%dx = gfld%igdtmpl(17) ! Di - i direction increment
404 map%dy = gfld%igdtmpl(18) ! N - # of parallels between pole and equator
405 map%lat1 = gfld%igdtmpl(12) ! La1 - lat of 1st grid point
406 map%lon1 = gfld%igdtmpl(13) ! Lo1 - lon of 1st grid point
407 write(tmp8,'(b8.8)') gfld%igdtmpl(14) ! resolution/component flag
408 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
409 map%r_earth = earth_radius (gfld%igdtmpl(1),
410 & gfld%igdtmpl(2),gfld%igdtmpl(3))
412 ! Scale dx/dy values to degrees, default range is 1e6.
413 if (map%dx.gt.10000) then
414 map%dx = map%dx/scale_factor
415 endif
416 if (map%dy.gt.10000) then
417 map%dy = (map%dy/scale_factor)*(-1)
418 endif
420 ! Scale lat/lon values to 0-180, default range is 1e6.
421 if (map%lat1.ge.scale_factor) then
422 map%lat1 = map%lat1/scale_factor
423 endif
424 if (map%lon1.ge.scale_factor) then
425 map%lon1 = map%lon1/scale_factor
426 endif
427 if ( debug_level .gt. 2 ) then
428 call mprintf(.true.,DEBUG,
429 & "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
430 & newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
431 & i1=nint(map%dy))
432 end if
434 else
435 call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
436 & newline=.true.,i1=gfld%igdtnum)
437 call mprintf(.true.,STDOUT,
438 & "see Code Table 3.1: Grid Definition Template Number",
439 & newline=.true.)
440 call mprintf(.true.,LOGFILE,
441 & "GRIB2 Unknown Projection: %i",
442 & newline=.true.,i1=gfld%igdtnum)
443 call mprintf(.true.,LOGFILE,
444 & "see Code Table 3.1: Grid Definition Template Number",
445 & newline=.true.)
446 endif
448 call mprintf(.true.,DEBUG,"G2 igrid = %i , dx = %f , dy = %
449 &f ", newline=.true., i1 = map%igrid, f1=map%dx, f2=map%dy)
451 if (icenter.eq.7) then
452 call ncep_grid_num (gfld%igdtnum)
453 endif
455 ! Deallocate arrays decoding GRIB2 record.
456 call gf_free(gfld)
457 endif
459 ! ----
461 ! Continue to unpack GRIB2 field.
462 do n=1,numfields ! e.g. U and V would =2, otherwise its usually =1
463 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
464 if (ierr.ne.0) then
465 write(*,*) ' ERROR extracting field gf_getfld = ',ierr
466 cycle
467 endif
469 ! The JMA GSM has two different grids in the same GRIB file, so we need
470 ! to process the map info for each field separately. If any other centers do
471 ! this, then processing will need to be added here, too.
473 if (icenter .eq. 34 .and. gfld%igdtnum.eq.0) then
474 map%nx = gfld%igdtmpl(8)
475 map%ny = gfld%igdtmpl(9)
476 map%dx = gfld%igdtmpl(17)
477 map%dy = gfld%igdtmpl(18)
478 ! Scale dx/dy values to degrees, default range is 1e6.
479 if (map%dx.gt.10000) then
480 map%dx = map%dx/scale_factor
481 endif
482 if (map%dy.gt.10000) then
483 map%dy = map%dy/scale_factor
484 endif
485 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
486 read(tmp8,'(1x,i1)') jscan
487 write(0,*) 'gfld%igdtmpl(19) = ',gfld%igdtmpl(19),
488 & ' jscan = ',jscan
489 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
490 map%dy = -1. * map%dy
491 endif
492 endif ! JMA spectral
494 ! ------------------------------------
495 ! Additional print information for developer.
496 if ( debug_level .GT. 1000 ) then
497 !MGD print *
498 !MGD print *,'G2 FIELD ',n
499 !MGD if (n==1) then
500 !MGD print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
501 !MGD print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
502 !MGD endif
503 !MGD if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
504 !MGD print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
505 !MGD endif
506 !MGD print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
507 !MGD & gfld%numoct_opt,gfld%interp_opt,
508 !MGD & gfld%igdtnum
509 !MGD print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
510 !MGD & (gfld%igdtmpl(j),j=1,gfld%igdtlen)
511 !MGD if ( gfld%num_opt .eq. 0 ) then
512 !MGD print *,'G2 NO Section 3 List Defining No. of Data Points.'
513 !MGD else
514 !MGD print *,'G2 Section 3 Optional List: ',
515 !MGD & (gfld%list_opt(j),j=1,gfld%num_opt)
516 !MGD endif
517 !MGD print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
518 !MGD & (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
520 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
521 & gfld%ipdtmpl(2))
522 !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
523 !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
524 !MGD print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
526 !MGD if ( gfld%num_coord .eq. 0 ) then
527 !MGD print *,'G2 NO Optional Vertical Coordinate List.'
528 !MGD else
529 !MGD print *,'G2 Section 4 Optional Coordinates: ',
530 !MGD & (gfld%coord_list(j),j=1,gfld%num_coord)
531 !MGD endif
532 if ( gfld%ibmap .ne. 255 ) then
533 call mprintf(.true.,DEBUG,
534 & 'G2 Num. of Data Points = %i with BIT-MAP %i',
535 & newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
536 else
537 call mprintf(.true.,DEBUG,
538 & 'G2 Num. of Data Points = %i NO BIT-MAP',
539 & newline=.true., i1=gfld%ndpts)
540 endif
541 !MGD print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
542 !MGD & (gfld%idrtmpl(j),j=1,gfld%idrtlen)
543 fldmax=gfld%fld(1)
544 fldmin=gfld%fld(1)
545 sum=gfld%fld(1)
546 do j=2,gfld%ndpts
547 if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
548 if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
549 sum=sum+gfld%fld(j)
550 enddo ! gfld%ndpts
552 !MGD print *,'G2 Data Values:'
553 call mprintf(.true.,DEBUG,'G2 MIN=%f AVE=%f MAX=%f',
554 & newline=.true., f1=fldmin, f2=sum/gfld%ndpts, f3=fldmax)
555 !do j=1,gfld%ndpts\20
556 ! write(*,*) j, gfld%fld(j)
557 !enddo
558 endif ! Additional Print information
559 ! ------------------------------------
561 ! do i = 1, maxvar
562 ! write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
563 ! enddo
564 !MGD if (debug_level .gt. 50) then
565 !MGD write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
566 !MGD & gfld%ipdtmpl(2),gfld%ipdtmpl(10)
567 !MGD endif
568 call mprintf(.true.,DEBUG,"G2 Searching the g2code array (Vta
569 &ble) for this grib field %i %i %i %i ", newline=.true.,
570 & i1 = gfld%discipline, i2 = gfld%ipdtmpl(1),
571 & i3 = gfld%ipdtmpl(2), i4 = gfld%ipdtmpl(10) )
573 ! Test this data record against list of desired variables
574 ! found in Vtable.
575 ! ----
576 MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
577 ! maxvar is defined in table.mod
579 if ( gfld%discipline .eq. g2code(1,i) .and. !Discipline
580 & gfld%ipdtmpl(1) .eq. g2code(2,i) .and. !Category
581 & gfld%ipdtmpl(2) .eq. g2code(3,i) .and. !Parameter
582 & gfld%ipdtmpl(10) .eq. g2code(4,i) .and. !Elevation
583 & gfld%ipdtnum .eq. g2code(5,i)) then !Template
585 call gf_free(gfld)
586 call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
587 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
588 & gfld%ipdtmpl(2))
590 !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
591 my_field=namvar(i)
593 !MGD if (debug_level .gt. 50) then
594 !MGD write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
595 !MGD write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
596 !MGD endif
599 ! need to match up soil levels with those requested.
600 ! For the Vtable levels, -88 = all levels, -99 = missing. The units
601 ! vary depending on the level code (e.g. 106 = cm, 103 = m).
602 if ( gfld%ipdtmpl(10) .eq. 106 ) then
603 TMP8LOOP: do j = 1, maxvar
604 if ((g2code(4,j) .eq. 106) .and.
605 & (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
606 & (gfld%ipdtmpl(12) .eq. level1(j)) .and.
607 & ((gfld%ipdtmpl(15) .eq. level2(j)) .or.
608 & (level2(j) .le. -88))) then
609 my_field = namvar(j)
610 exit TMP8LOOP
611 endif
612 enddo TMP8LOOP
613 if (j .gt. maxvar ) then
614 write(6,'(a,i6,a,i6,a)') 'Subsoil level ',
615 & gfld%ipdtmpl(12),' to ',gfld%ipdtmpl(15),
616 & ' in the GRIB2 file, was not found in the Vtable'
617 cycle MATCH_LOOP
618 endif
619 !MGD if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
620 endif
622 ! Level (eg. 10000 mb)
623 if(gfld%ipdtmpl(10).eq.100) then
624 ! Pressure level (range from 1000mb to 0mb)
625 level=gfld%ipdtmpl(12) *
626 & (10. ** (-1. * gfld%ipdtmpl(11)))
627 elseif(gfld%ipdtmpl(10).eq.105) then
628 ! Hybrid level (range from 1 to N)
629 level=gfld%ipdtmpl(12)
630 elseif(gfld%ipdtmpl(10).eq.104) then
631 ! Sigma level (range from 10000 to 0)
632 level=gfld%ipdtmpl(12)
633 elseif(gfld%ipdtmpl(10).eq.101) then
634 ! MSL
635 level=201300.
636 elseif(gfld%ipdtmpl(10).eq.106.or.
637 & gfld%ipdtmpl(10).eq.1) then
638 ! Misc near ground/surface levels
639 level=200100.
640 elseif(gfld%ipdtmpl(10).eq.6) then
641 ! Level of Max wind
642 level=6. ! .06 mb should never be used by anyone
643 elseif(gfld%ipdtmpl(10).eq.7) then
644 ! Tropopause
645 level=7. ! .07 mb should never be used by anyone
646 else
647 ! Misc near ground/surface levels
648 level=200100.
649 endif
650 iplvl = int(level)
652 ! Store the unpacked 2D slab from the GRIB2 record
653 allocate(hold_array(gfld%ngrdpts))
654 do j=1,gfld%ngrdpts
655 hold_array(j)=gfld%fld(j)
656 enddo
658 ! Some grids need to be reordered. Until we get an example, this is
659 ! a placeholder
660 ! call reorder_it (hold_array, map%nx, map%ny, map%dx,
661 ! & map%dy, iorder)
663 ! When we have reached this point, we have a data array ARRAY
664 ! which has some data we want to save, with field name FIELD
665 ! at pressure level LEVEL (Pa). Dimensions of this data are
666 ! map%nx and map%ny. Put this data into storage.
668 !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
669 !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
670 call put_storage(iplvl,my_field,
671 & reshape(hold_array(1:map%nx*map%ny),
672 & (/map%nx, map%ny/)), map%nx,map%ny)
673 deallocate(hold_array)
675 ! If Specific Humidity is present on hybrid levels AND
676 ! upper-air RH is missing, see if we can compute RH from
677 ! Specific Humidity.
678 if (.not. is_there(iplvl, 'RH') .and.
679 & is_there(iplvl, 'SH') .and.
680 & is_there(iplvl, 'TT') .and.
681 & is_there(iplvl, 'P')) then
682 call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
683 !call llstor_remove(iplvl, 'SH') !We are done with SH
684 endif
686 ith=ith+1
687 exit MATCH_LOOP
689 endif ! Selected param.
692 enddo MATCH_LOOP
694 ! Deallocate arrays decoding GRIB2 record.
695 call gf_free(gfld)
697 enddo ! 1,numfields
700 enddo VERSION ! skgb
703 if ( debug_level .gt. 100 ) then
704 call mprintf (.true.,DEBUG,
705 & "G2 total number of fields found = %i ",newline=.true.,i1=itot)
706 end if
708 CALL BACLOSE(junit,IOS)
710 nullify(gfld%local) ! must be nullified before opening next file
711 ireaderr=1
712 else
713 call mprintf (.true.,DEBUG,"open status failed because %i ",
714 & newline=.true., i1=ios)
715 hdate = '9999-99-99_99:99:99'
716 ireaderr=2
717 endif ! ireaderr check
719 END subroutine rd_grib2
721 !*****************************************************************************!
722 ! Subroutine edition_num !
723 ! !
724 ! Purpose: !
725 ! Read one record from the input GRIB2 file. Based on the information in !
726 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
727 ! the GRIB2 record is one to process or to skip. If the field is one we !
728 ! want to keep, extract the data from the GRIB2 record, and pass the data !
729 ! back to the calling routine. !
730 ! !
731 ! Argument list: !
732 ! Input: !
733 ! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
734 ! unit number, since we do not do Fortran I/O for the GRIB2 !
735 ! files. Nor is it a UNIX File Descriptor returned from a C !
736 ! OPEN statement. It is really just an array index to the !
737 ! array (IUARR) where the UNIX File Descriptor values are !
738 ! stored. !
739 ! GRIB2FILE: File name to open, if it is not already open. !
740 ! !
741 ! Output: !
742 ! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
743 ! IERR : Error flag: 0 - no error on read from GRIB2 file. !
744 ! 1 - Hit the end of the GRIB2 file. !
745 ! 2 - The file GRIBFLNM we tried to open does !
746 ! not exist. !
747 ! Author: Paula McCaslin !
748 ! NOAA/FSL !
749 ! Sept 2004 !
750 !*****************************************************************************!
752 SUBROUTINE edition_num(junit, gribflnm,
753 & grib_edition, ireaderr)
755 use grib_mod
756 use params
757 use module_debug
759 parameter(msk1=32000,msk2=4000)
760 character(len=1),allocatable,dimension(:) :: cgrib
761 integer :: listsec0(3)
762 integer :: listsec1(13)
763 character(len=*) :: gribflnm
764 integer :: lskip, lgrib
765 integer :: junit
766 integer :: grib_edition
767 integer :: i, j, ireaderr
768 integer :: currlen
770 character(len=4) :: ctemp
771 character(len=4),parameter :: grib='GRIB',c7777='7777'
773 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
774 C SET ARGUMENTS
776 itot=0
777 icount=0
778 iseek=0
779 lskip=0
780 lgrib=0
781 currlen=0
783 !/* IOS Return Codes from BACIO: */
784 !/* 0 All was well */
785 !/* -1 Tried to open read only _and_ write only */
786 !/* -2 Tried to read and write in the same call */
787 !/* -3 Internal failure in name processing */
788 !/* -4 Failure in opening file */
789 !/* -5 Tried to read on a write-only file */
790 !/* -6 Failed in read to find the 'start' location */
791 !/* -7 Tried to write to a read only file */
792 !/* -8 Failed in write to find the 'start' location */
793 !/* -9 Error in close */
794 !/* -10 Read or wrote fewer data than requested */
796 !if ireaderr =1 we have hit the end of a file.
797 !if ireaderr =2 we have hit the end of all the files.
798 !if ireaderr =3 beginning characters 'GRIB' not found
800 ! Open a byte-addressable file.
801 CALL BAOPENR(junit,gribflnm,IOS)
802 if (ios.eq.0) then
804 ! Search opend file for the next GRIB2 messege (record).
805 call skgb(junit,iseek,msk1,lskip,lgrib)
807 ! Check for EOF, or problem
808 call mprintf((lgrib.eq.0),ERROR,
809 & "Grib2 file or date problem, stopping in edition_num.",
810 & newline=.true.)
812 ! Check size, if needed allocate more memory.
813 if (lgrib.gt.currlen) then
814 if (allocated(cgrib)) deallocate(cgrib)
815 allocate(cgrib(lgrib),stat=is)
816 currlen=lgrib
817 endif
819 ! Read a given number of bytes from unblocked file.
820 call baread(junit,lskip,lgrib,lengrib,cgrib)
822 ! Check for beginning of GRIB message in the first 100 bytes
823 istart=0
824 do j=1,100
825 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
826 if (ctemp.eq.grib ) then
827 istart=j
828 exit
829 endif
830 enddo
831 if (istart.eq.0) then
832 ireaderr=3
833 print*, "The beginning 4 characters >GRIB< were not found."
834 endif
836 ! Unpack Section 0 - Indicator Section to extract GRIB edition field
837 iofst=8*(istart+5)
838 call gbyte(cgrib,discipline,iofst,8) ! Discipline
839 iofst=iofst+8
840 call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
842 print *, 'ungrib - grib edition num', grib_edition
843 CALL BACLOSE(junit,IOS)
844 ireaderr=1
845 else if (ios .eq. -4) then
846 call mprintf(.true.,ERROR,
847 & "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
848 else
849 print *,'edition_num: open status failed because',ios,gribflnm
850 ireaderr=2
851 endif ! ireaderr check
853 END subroutine edition_num
855 !*****************************************************************************!
857 SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
858 ! Compute relative humidity from specific humidity in the upper air.
859 use storage_module
860 implicit none
861 integer :: ix, jx
862 integer :: iiplvl
863 real :: lat1, lon1, dx, dy
864 real, dimension(ix,jx) :: T, P, RH, Q
866 real, parameter :: svp1=611.2
867 real, parameter :: svp2=17.67
868 real, parameter :: svp3=29.65
869 real, parameter :: svpt0=273.15
870 real, parameter :: eps = 0.622
872 real startlat, startlon, deltalat, deltalon
874 call get_storage(iiplvl, 'P', P, ix, jx)
875 call get_storage(iiplvl, 'TT', T, ix, jx)
876 call get_storage(iiplvl, 'SH', Q, ix, jx)
878 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
880 call put_storage(iiplvl, 'RH', rh, ix, jx)
882 end subroutine g2_compute_rh_spechumd_upa
884 !*****************************************************************************!
886 subroutine ncep_grid_num (pnum)
887 !
888 ! Grib2 doesn't have a grid-number entry, so we have to figure it out
889 !
890 use gridinfo ! Included to define map%
891 integer :: pnum
892 real, parameter :: eps = .01
893 character (len=8) :: tmp8
895 ! write(6,*) 'begin ncep_grid_num'
896 ! write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
897 tmp8 = ' '
898 if (pnum .eq. 30) then
899 if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
900 write(tmp8,'("GRID 218")')
901 else if (abs(map%dx - 40.63525) .lt. eps
902 & .and. map%nx .eq. 185) then
903 write(tmp8,'("GRID 212")')
904 else if (abs(map%dx - 40.63525) .lt. eps
905 & .and. map%nx .eq. 151) then
906 write(tmp8,'("GRID 236")')
907 else if (abs(map%dx - 81.2705) .lt. eps
908 & .and. map%nx .eq. 93) then
909 write(tmp8,'("GRID 211")')
910 else if (abs (map%dx - 32.46341) .lt. eps
911 & .and. map%nx .eq. 349) then
912 write(tmp8,'("GRID 221")')
913 else if (abs(map%dx - 20.317625) .lt. eps
914 & .and. map%nx .eq. 301) then
915 write(tmp8,'("GRID 252")')
916 endif
917 else if (pnum .eq. 20) then
918 if (abs(map%dx - 15.0) .lt. eps) then
919 write(tmp8,'("GRID 88")')
920 endif
921 else if (pnum .eq. 0) then
922 if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
923 write(tmp8,'("GRID 3")')
924 else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
925 write(tmp8,'("GRID 4")')
926 endif
927 endif
928 map%source(25:32) = tmp8
929 ! write(6,*) 'map%source = ',map%source
930 end subroutine ncep_grid_num
931 !*****************************************************************************!
933 function earth_radius (icode, iscale, irad_m)
934 ! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
935 use module_debug
936 real :: earth_radius
937 integer :: icode
938 integer :: iscale, irad_m
939 if ( icode .eq. 0 ) then
940 earth_radius = 6367470. * .001
941 else if ( icode .eq. 1) then
942 earth_radius = 0.001 * float(irad_m) / 10**iscale
943 else if ( icode .eq. 6 ) then
944 earth_radius = 6371229. * .001
945 else if ( icode .eq. 8 ) then
946 earth_radius = 6371200. * .001
947 else
948 call mprintf(.true.,ERROR,
949 & "unknown earth radius for code %i",newline=.true.,i1=icode)
950 endif
951 end function earth_radius