| blob | 2b6d8d80fb687356f2883d44aa98908fbad7df18 |
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 ! pmin : Minimum pressure level (Pa) to process. !
22 ! !
23 ! Output: !
24 ! !
25 ! hdate : The (up to)19-character date of the field to process. !
26 ! grib_edition : Version of the gribfile (1 or 2) !
27 ! ireaderr : Error flag: 0 - no error on read from GRIB2 file. !
28 ! 1 - Hit the end of the GRIB2 file. !
29 ! 2 - The file GRIBFLNM we tried to open does !
30 ! not exist. !
31 ! !
32 ! !
33 ! Author: Paula McCaslin, NOAA/FSL, Sept 2004 !
34 ! Code is based on code developed by Steve Gilbert NCEP & Kevin Manning NCAR !
35 ! Adapted for WPS: Jim Bresch, NCAR/MMM. Sept 2006 !
36 !*****************************************************************************!
38 SUBROUTINE rd_grib2(junit, gribflnm, hdate,
39 & grib_edition, ireaderr, debug_level, pmin)
41 use grib_mod
42 use params
43 use table ! Included to define g2code
44 use gridinfo ! Included to define map%
45 use storage_module ! Included sub put_storage
46 use module_debug
48 real, allocatable, dimension(:) :: hold_array
49 parameter(msk1=32000,msk2=4000)
50 character(len=1),allocatable,dimension(:) :: cgrib
51 integer :: listsec0(3)
52 integer :: listsec1(13)
53 integer year, month, day, hour, minute, second, fcst
54 character(len=*) :: gribflnm
55 character(len=*) :: hdate
56 character(len=8) :: pabbrev
57 character(len=20) :: labbrev
58 character(len=80) :: tabbrev
59 integer :: lskip, lgrib
60 integer :: junit, itot, icount, iseek
61 integer :: grib_edition
62 integer :: i, j, ireaderr, ith , debug_level
63 integer :: currlen
64 logical :: unpack, expand
65 type(gribfield) :: gfld
66 ! For subroutine put_storage
67 real :: level
68 real :: scale_factor
69 integer :: iplvl
70 character (len=9) :: my_field
71 character (len=8) :: tmp8
72 ! For subroutine output
73 integer , parameter :: maxlvl = 150
74 real , dimension(maxlvl) :: plvl
75 integer :: nlvl
76 integer , dimension(maxlvl) :: level_array
77 real :: glevel1, glevel2
78 logical :: first = .true.
79 real :: pmin
81 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 C SET ARGUMENTS
84 unpack=.true.
85 expand=.true.
86 hdate = '0000-00-00_00:00:00'
87 ierr=0
88 itot=0
89 icount=0
90 iseek=0
91 lskip=0
92 lgrib=0
93 currlen=0
94 ith=1
95 scale_factor = 1e6
96 call mprintf(.true.,DEBUG,"Begin rd_grib2", newline=.true.)
98 !/* IOS Return Codes from BACIO: */
99 !/* 0 All was well */
100 !/* -1 Tried to open read only _and_ write only */
101 !/* -2 Tried to read and write in the same call */
102 !/* -3 Internal failure in name processing */
103 !/* -4 Failure in opening file */
104 !/* -5 Tried to read on a write-only file */
105 !/* -6 Failed in read to find the 'start' location */
106 !/* -7 Tried to write to a read only file */
107 !/* -8 Failed in write to find the 'start' location */
108 !/* -9 Error in close */
109 !/* -10 Read or wrote fewer data than requested */
111 !if ireaderr =1 we have hit the end of a file.
112 !if ireaderr =2 we have hit the end of all the files.
115 ! Open a byte-addressable file.
116 CALL BAOPENR(junit,gribflnm,IOS)
117 first = .true.
118 if (ios.eq.0) then
119 VERSION: do
121 ! Search opend file for the next GRIB2 messege (record).
122 call skgb(junit,iseek,msk1,lskip,lgrib)
124 ! Check for EOF, or problem
125 if (lgrib.eq.0) then
126 exit
127 endif
129 ! Check size, if needed allocate more memory.
130 if (lgrib.gt.currlen) then
131 if (allocated(cgrib)) deallocate(cgrib)
132 allocate(cgrib(lgrib),stat=is)
133 !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
134 currlen=lgrib
135 endif
137 ! Read a given number of bytes from unblocked file.
138 call baread(junit,lskip,lgrib,lengrib,cgrib)
140 call mprintf ((lgrib.ne.lengrib),ERROR,
141 & "rd_grib2: IO Error. %i .ne. %i ", newline=.true.,
142 & i1=lgrib,i2=lengrib)
144 iseek=lskip+lgrib
145 icount=icount+1
147 call mprintf (.true.,DEBUG,
148 & "G2 GRIB MESSAGE %i starts at %i ", newline=.true.,
149 & i1=icount,i2=lskip+1)
151 ! Unpack GRIB2 field
152 call gb_info(cgrib,lengrib,listsec0,listsec1,
153 & numfields,numlocal,maxlocal,ierr)
154 call mprintf((ierr.ne.0),ERROR,
155 & " ERROR querying GRIB2 message = %i",newline=.true.,i1=ierr)
156 itot=itot+numfields
158 grib_edition=listsec0(2)
159 if (grib_edition.ne.2) then
160 exit VERSION
161 endif
163 ! Additional print statments for developer.
164 !MGD if ( debug_level .GT. 100 ) then
165 !MGD print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
166 !MGD print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
167 !MGD print *,'G2 Contains ',numlocal,' Local Sections ',
168 !MGD & ' and ',numfields,' data fields.'
169 !MGD endif
172 ! ----
173 ! Once per file fill in date, model and projection values.
175 if (first) then
176 first = .false.
178 ! Build the 19-character date string, based on GRIB2 header date
179 ! and time information, including forecast time information:
181 n=1
182 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
185 year =gfld%idsect(6) !(FOUR-DIGIT) YEAR OF THE DATA
186 month =gfld%idsect(7) ! MONTH OF THE DATA
187 day =gfld%idsect(8) ! DAY OF THE DATA
188 hour =gfld%idsect(9) ! HOUR OF THE DATA
189 minute=gfld%idsect(10) ! MINUTE OF THE DATA
190 second=gfld%idsect(11) ! SECOND OF THE DATA
192 fcst = 0
194 ! Extract forecast time. Assume the first field's valid time is true for all fields.
195 ! This doesn't have to be true, but ungrib is designed to decode one time-level at
196 ! a time.
198 if ( gfld%ipdtnum .ne. 8 ) then
199 if ( gfld%ipdtmpl(8) .eq. 1 ) then ! time units are hours
200 fcst = gfld%ipdtmpl(9)
201 else if ( gfld%ipdtmpl(8) .eq. 0 ) then ! minutes
202 fcst = gfld%ipdtmpl(9) / 60.
203 else if ( gfld%ipdtmpl(8) .eq. 2 ) then ! days
204 fcst = gfld%ipdtmpl(9) * 24.
205 else
206 call mprintf(.true.,ERROR,
207 & "Time unit in ipdtmpl(8), %i is not suported",
208 & newline=.true.,i1=gfld%ipdtmpl(8))
209 endif
210 else
211 ! pdt 4.8 data are time-averaged, accumulated, or min/max fields with the
212 ! ending (valid) time provided.
213 year =gfld%ipdtmpl(16)
214 month =gfld%ipdtmpl(17)
215 day =gfld%ipdtmpl(18)
216 hour =gfld%ipdtmpl(19)
217 minute=gfld%ipdtmpl(20)
218 second=gfld%ipdtmpl(21)
219 fcst = 0.
220 endif
222 if ( gfld%idsect(5) .eq. 2 ) fcst = 0.
223 ! Compute valid time.
225 !print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
226 !print *, 'hhmm ',gfld%idsect(9),gfld%idsect(10)
228 call build_hdate(hdate,year,month,day,hour,minute,second)
229 call mprintf(.true.,DEBUG,"G2 hdate = %s ", newline=.true.,
230 & s1=hdate)
231 call geth_newdate(hdate,hdate,3600*fcst)
232 call mprintf(.true.,DEBUG,"G2 hdate (fcst?) = %s ",
233 & newline=.true., s1=hdate)
235 !--
237 ! Indicator of the source (center) of the data.
238 icenter = gfld%idsect(1)
240 ! Indicator of model (or whatever) which generated the data.
241 iprocess = gfld%ipdtmpl(5)
244 if (icenter.eq.7) then
245 if (iprocess.eq.81) then
246 map%source = 'NCEP GFS Analysis'
247 elseif (iprocess.eq.82) then
248 map%source = 'NCEP GFS GDAS/FNL'
249 elseif (iprocess.eq.83) then
250 map%source = 'NCEP HRRR Model'
251 elseif (iprocess.eq.84) then
252 map%source = 'NCEP MESO NAM Model'
253 elseif (iprocess.eq.89) then
254 map%source = 'NCEP NMM '
255 elseif (iprocess.eq.96) then
256 map%source = 'NCEP GFS Model'
257 elseif (iprocess.eq.86 .or. iprocess.eq.100) then
258 map%source = 'NCEP RUC Model' ! 60 km
259 elseif (iprocess.eq.101) then
260 map%source = 'NCEP RUC Model' ! 40 km
261 elseif (iprocess.eq.105) then
262 if (year .gt. 2011) then
263 map%source = 'NCEP RAP Model'
264 else
265 map%source = 'NCEP RUC Model' ! 20 km
266 endif
267 elseif (iprocess.eq.107) then
268 map%source = 'NCEP GEFS'
269 elseif (iprocess.eq.109) then
270 map%source = 'NCEP RTMA'
271 elseif (iprocess.eq.140) then
272 map%source = 'NCEP NARR'
273 elseif (iprocess.eq.44) then
274 map%source = 'NCEP SST Analysis'
275 elseif (iprocess.eq.70) then
276 map%source = 'GFDL Hurricane Model'
277 elseif (iprocess.eq.80) then
278 map%source = 'NCEP GFS Ensemble'
279 elseif (iprocess.eq.107) then ! renumbered as of 23 Feb 2010
280 map%source = 'NCEP GFS Ensemble'
281 elseif (iprocess.eq.111) then
282 map%source = 'NCEP NMMB Model'
283 elseif (iprocess.eq.112) then
284 map%source = 'NCEP WRF-NMM Model'
285 elseif (iprocess.eq.116) then
286 map%source = 'NCEP WRF-ARW Model'
287 elseif (iprocess.eq.129) then
288 map%source = 'NCEP GODAS'
289 elseif (iprocess.eq.197) then
290 map%source = 'NCEP CDAS CFSV2'
291 elseif (iprocess.eq.25) then
292 map%source = 'NCEP SNOW COVER ANALYSIS'
293 else
294 map%source = 'unknown model from NCEP'
295 call mprintf(.true.,STDOUT,
296 & "unknown model from NCEP %i ",newline=.true.,
297 & i1=iprocess)
298 call mprintf(.true.,LOGFILE,
299 & "unknown model from NCEP %i ",newline=.true.,
300 & i1=iprocess)
301 end if
302 else if (icenter .eq. 57) then
303 if (iprocess .eq. 87) then
304 map%source = 'AFWA AGRMET'
305 else
306 map%source = 'AFWA'
307 endif
308 else if ( icenter .eq. 58 ) then
309 map%source = 'US Navy FNOC'
310 else if (icenter .eq. 59) then
311 if (iprocess .eq. 125) then
312 map%source = 'NOAA GSD Rapid Refresh Model'
313 else if (iprocess .eq. 83) then
314 map%source = 'NOAA GSD HRRR Model'
315 else if (iprocess .eq. 105) then
316 map%source = 'NOAA GSD'
317 else
318 print *,'Unknown GSD source'
319 stop
320 endif
321 else if (icenter .eq. 60) then
322 map%source = 'NCAR'
323 else if (icenter .eq. 98) then
324 map%source = 'ECMWF'
325 else if (icenter .eq. 34) then
326 map%source = 'JMA'
327 else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
328 map%source = 'UKMO'
329 else
330 map%source = 'unknown model and orig center'
331 end if
332 if (icenter.eq.7) then
333 if (iprocess.eq.81 .or. iprocess.eq.82 .or.
334 & iprocess.eq.96) then
335 ! pmin should not be 15 or 40 hPa for gfs files. Use the next highest level.
336 if (pmin .eq. 1500.) then
337 pmin = 1000.
338 else if (pmin .eq. 4000.) then
339 pmin = 3000.
340 endif
341 endif
342 endif
344 call mprintf(.true.,DEBUG,"G2 source is = %s ",
345 & newline=.true., s1=map%source)
347 !--
349 ! Store information about the grid containing the data.
350 ! This stuff gets stored in the MAP variable, as defined in
351 ! module GRIDINFO.
353 map%startloc = 'SWCORNER'
354 map%grid_wind = .true.
356 if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
357 map%igrid = 0
358 map%nx = gfld%igdtmpl(8)
359 map%ny = gfld%igdtmpl(9)
360 map%dx = gfld%igdtmpl(17)
361 map%dy = gfld%igdtmpl(18)
362 map%lat1 = gfld%igdtmpl(12)
363 map%lon1 = gfld%igdtmpl(13)
364 write(tmp8,'(b8.8)') gfld%igdtmpl(14)
365 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
366 map%r_earth = earth_radius (gfld%igdtmpl(1),
367 & gfld%igdtmpl(2),gfld%igdtmpl(3))
369 ! Fix for NCEP 1/12 degree grids (e.g. rtgsst)
370 if (icenter .eq. 7 .and. map%dx .eq. 83000. .and. map%nx
371 & .eq. 4320) then
372 map%lat1 = 89958333.
373 map%lon1 = 41667.
374 map%dx = 83333.333 * sign(1.0,map%dx)
375 map%dy = 83333.333 * sign(1.0,map%dy)
376 endif
378 if ((gfld%igdtmpl(10) .eq. 0).OR.
379 & (gfld%igdtmpl(10) .eq. 255)) THEN
380 ! Scale lat/lon values to 0-180, default range is 1e6.
381 map%lat1 = map%lat1/scale_factor
382 map%lon1 = map%lon1/scale_factor
383 ! Scale dx/dy values to degrees, default range is 1e6.
384 map%dx = map%dx/scale_factor
385 map%dy = map%dy/scale_factor
386 else
387 ! Basic angle and subdivisions are non-zero (not tested)
388 map%lat1 = map%lat1 *
389 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
390 map%lon1 = map%lon1 *
391 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
392 map%dx = map%dx *
393 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
394 map%dy = map%dy *
395 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
396 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
397 &has not been tested, continuing anyway")
398 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
399 & has not been tested, continuing anyway")
400 endif
403 ! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
404 ! In WPS, the sign of dy indicates the direction of the scan.
405 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
406 read(tmp8,'(1x,i1)') jscan
407 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
408 map%dy = -1. * map%dy
409 endif
410 ! if ( map%lat1 .gt. gfld%igdtmpl(15) .and.
411 ! & map%dy .gt. 0. ) then
412 ! map%dy = -1. * map%dy
413 ! write(6,*) 'Resetting map%dy for iprocess = ',iprocess
414 ! endif
416 elseif (gfld%igdtnum.eq.10) then ! Mercator Grid.
417 map%igrid = 1
418 map%nx = gfld%igdtmpl(8)
419 map%ny = gfld%igdtmpl(9)
420 map%lov = 0.
421 map%truelat1 = gfld%igdtmpl(13) / scale_factor
422 map%truelat2 = 0.
423 map%dx = gfld%igdtmpl(18) / scale_factor
424 map%dy = gfld%igdtmpl(19) / scale_factor
425 map%lat1 = gfld%igdtmpl(10) / scale_factor
426 map%lon1 = gfld%igdtmpl(11) / scale_factor
427 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
428 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
429 map%r_earth = earth_radius (gfld%igdtmpl(1),
430 & gfld%igdtmpl(2),gfld%igdtmpl(3))
432 elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
433 map%igrid = 5
434 map%nx = gfld%igdtmpl(8)
435 map%ny = gfld%igdtmpl(9)
436 map%lov = gfld%igdtmpl(14) / scale_factor
437 map%truelat1 = 60.
438 map%truelat2 = 91.
439 map%dx = gfld%igdtmpl(15) / scale_factor
440 map%dy = gfld%igdtmpl(16) / scale_factor
441 map%lat1 = gfld%igdtmpl(10) / scale_factor
442 map%lon1 = gfld%igdtmpl(11) / scale_factor
443 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
444 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
445 map%r_earth = earth_radius (gfld%igdtmpl(1),
446 & gfld%igdtmpl(2),gfld%igdtmpl(3))
448 elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
449 map%igrid = 3
450 map%nx = gfld%igdtmpl(8)
451 map%ny = gfld%igdtmpl(9)
452 map%lov = gfld%igdtmpl(14) / scale_factor
453 map%truelat1 = gfld%igdtmpl(19) / scale_factor
454 map%truelat2 = gfld%igdtmpl(20) / scale_factor
455 map%dx = gfld%igdtmpl(15) / scale_factor
456 map%dy = gfld%igdtmpl(16) / scale_factor
457 map%lat1 = gfld%igdtmpl(10) / scale_factor
458 map%lon1 = gfld%igdtmpl(11) / scale_factor
459 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
460 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
461 map%r_earth = earth_radius (gfld%igdtmpl(1),
462 & gfld%igdtmpl(2),gfld%igdtmpl(3))
464 elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
465 map%igrid = 4
466 map%nx = gfld%igdtmpl(8) ! Ni - # of points along a parallel
467 map%ny = gfld%igdtmpl(9) ! Nj - # of points along meridian
468 map%dx = gfld%igdtmpl(17) ! Di - i direction increment
469 map%dy = gfld%igdtmpl(18) ! N - # of parallels between pole and equator
470 map%lat1 = gfld%igdtmpl(12) ! La1 - lat of 1st grid point
471 map%lon1 = gfld%igdtmpl(13) ! Lo1 - lon of 1st grid point
472 write(tmp8,'(b8.8)') gfld%igdtmpl(14) ! resolution/component flag
473 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
474 map%r_earth = earth_radius (gfld%igdtmpl(1),
475 & gfld%igdtmpl(2),gfld%igdtmpl(3))
477 ! Scale dx/dy values to degrees, default range is 1e6.
478 if (map%dx.gt.10000) then
479 map%dx = map%dx/scale_factor
480 endif
481 if (map%dy.gt.10000) then
482 map%dy = (map%dy/scale_factor)*(-1)
483 endif
485 ! Fix for zonal shift in CFSR data, following a similar fix
486 ! for global lat-lon data in rd_grib1.F
487 if ( ABS(map%nx * map%dx - 360.0) < 1.0 ) then
488 if (ABS(map%dx - (360./real(map%nx))) > 0.00001) then
489 write(0,*) 'CFSR fix: recomputing delta-longitude'
490 map%dx = 360./real(map%nx)
491 endif
492 endif
494 ! Scale lat/lon values to 0-180, default range is 1e6.
495 if (abs(map%lat1).ge.scale_factor) then
496 map%lat1 = map%lat1/scale_factor
497 endif
498 if (map%lon1.ge.scale_factor) then
499 map%lon1 = map%lon1/scale_factor
500 endif
501 if ( debug_level .gt. 2 ) then
502 call mprintf(.true.,DEBUG,
503 & "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
504 & newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
505 & i1=nint(map%dy))
506 end if
507 elseif (gfld%igdtnum.eq.32769) then ! Arakawa Non-E Staggered grid.
508 map%igrid = 6
509 map%nx = gfld%igdtmpl(8)
510 map%ny = gfld%igdtmpl(9)
511 map%dx = gfld%igdtmpl(17)
512 map%dy = gfld%igdtmpl(18)
513 map%lat1 = gfld%igdtmpl(12)
514 map%lon1 = gfld%igdtmpl(13)
515 map%lat0 = gfld%igdtmpl(15)
516 map%lon0 = gfld%igdtmpl(16)
517 map%r_earth = earth_radius (gfld%igdtmpl(1),
518 & gfld%igdtmpl(2),gfld%igdtmpl(3))
519 if ((gfld%igdtmpl(10) .eq. 0).OR.
520 & (gfld%igdtmpl(10) .eq. 255)) THEN
521 map%lat1 = map%lat1/scale_factor
522 map%lon1 = map%lon1/scale_factor
523 map%lat0 = map%lat0/scale_factor
524 map%lon0 = map%lon0/scale_factor
525 map%dx = map%dx/scale_factor/1.e3
526 map%dy = map%dy/scale_factor/1.e3
527 else
528 ! Basic angle and subdivisions are non-zero (not tested)
529 map%lat1 = map%lat1 *
530 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
531 map%lon1 = map%lon1 *
532 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
533 map%dx = map%dx *
534 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
535 map%dy = map%dy *
536 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
537 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
538 &has not been tested, continuing anyway")
539 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
540 & has not been tested, continuing anyway")
541 endif
543 else
544 call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
545 & newline=.true.,i1=gfld%igdtnum)
546 call mprintf(.true.,STDOUT,
547 & "ungrib understands projections 0, 20, 30, and 40",
548 & newline=.true.)
549 call mprintf(.true.,LOGFILE,
550 & "GRIB2 Unknown Projection: %i",
551 & newline=.true.,i1=gfld%igdtnum)
552 call mprintf(.true.,LOGFILE,
553 & "ungrib understands projections 0, 10, 20, 30, and 40",
554 & newline=.true.)
555 ! If the projection is not known, then it can't be processed by metgrid/plotfmt
556 stop 'Stop in rd_grib2'
557 endif
559 call mprintf(.true.,DEBUG,"G2 igrid = %i , dx = %f , dy = %
560 &f ", newline=.true., i1 = map%igrid, f1=map%dx, f2=map%dy)
562 if (icenter.eq.7) then
563 call ncep_grid_num (gfld%igdtnum)
564 endif
566 ! Deallocate arrays decoding GRIB2 record.
567 call gf_free(gfld)
569 endif ! "first" if-block
571 ! ----
573 ! Continue to unpack GRIB2 field.
574 NUM_FIELDS: do n = 1, numfields
575 ! e.g. U and V would =2, otherwise its usually =1
576 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
577 if (ierr.ne.0) then
578 write(*,*) ' ERROR extracting field gf_getfld = ',ierr
579 cycle
580 endif
582 ! The JMA GSM has two different grids in the same GRIB file, so we need
583 ! to process the map info for each field separately. If any other centers do
584 ! this, then processing will need to be added here, too.
586 if (icenter .eq. 34 .and. gfld%igdtnum.eq.0) then
587 map%nx = gfld%igdtmpl(8)
588 map%ny = gfld%igdtmpl(9)
589 map%dx = gfld%igdtmpl(17)
590 map%dy = gfld%igdtmpl(18)
591 ! Scale dx/dy values to degrees, default range is 1e6.
592 if (map%dx.gt.10000) then
593 map%dx = map%dx/scale_factor
594 endif
595 if (map%dy.gt.10000) then
596 map%dy = map%dy/scale_factor
597 endif
598 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
599 read(tmp8,'(1x,i1)') jscan
600 write(0,*) 'gfld%igdtmpl(19) = ',gfld%igdtmpl(19),
601 & ' jscan = ',jscan
602 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
603 map%dy = -1. * map%dy
604 endif
605 endif ! JMA spectral
607 ! ------------------------------------
608 ! Additional print information for developer.
609 if ( debug_level .GT. 1000 ) then
610 !MGD print *
611 !MGD print *,'G2 FIELD ',n
612 !MGD if (n==1) then
613 !MGD print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
614 !MGD print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
615 !MGD endif
616 !MGD if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
617 !MGD print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
618 !MGD endif
619 !MGD print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
620 !MGD & gfld%numoct_opt,gfld%interp_opt,
621 !MGD & gfld%igdtnum
622 !MGD print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
623 !MGD & (gfld%igdtmpl(j),j=1,gfld%igdtlen)
624 !MGD if ( gfld%num_opt .eq. 0 ) then
625 !MGD print *,'G2 NO Section 3 List Defining No. of Data Points.'
626 !MGD else
627 !MGD print *,'G2 Section 3 Optional List: ',
628 !MGD & (gfld%list_opt(j),j=1,gfld%num_opt)
629 !MGD endif
630 !MGD print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
631 !MGD & (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
633 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
634 & gfld%ipdtmpl(2))
635 !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
636 !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
637 !MGD print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
639 !MGD if ( gfld%num_coord .eq. 0 ) then
640 !MGD print *,'G2 NO Optional Vertical Coordinate List.'
641 !MGD else
642 !MGD print *,'G2 Section 4 Optional Coordinates: ',
643 !MGD & (gfld%coord_list(j),j=1,gfld%num_coord)
644 !MGD endif
645 if ( gfld%ibmap .ne. 255 ) then
646 call mprintf(.true.,DEBUG,
647 & 'G2 Num. of Data Points = %i with BIT-MAP %i',
648 & newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
649 else
650 call mprintf(.true.,DEBUG,
651 & 'G2 Num. of Data Points = %i NO BIT-MAP',
652 & newline=.true., i1=gfld%ndpts)
653 endif
654 !MGD print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
655 !MGD & (gfld%idrtmpl(j),j=1,gfld%idrtlen)
656 endif ! Additional Print information
657 ! ------------------------------------
659 ! do i = 1, maxvar
660 ! write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
661 ! enddo
662 !MGD if (debug_level .gt. 50) then
663 !MGD write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
664 !MGD & gfld%ipdtmpl(2),gfld%ipdtmpl(10)
665 !MGD endif
666 call mprintf(.true.,DEBUG,"G2 Searching the g2code array (Vta
667 &ble) for this grib field %i %i %i %i %i %i ", newline=.true.,
668 & i1 = gfld%discipline, i2 = gfld%ipdtmpl(1),
669 & i3 = gfld%ipdtmpl(2), i4 = gfld%ipdtmpl(10),
670 & i5 = gfld%ipdtmpl(12), i6 = gfld%ipdtnum )
673 ! Test this data record against list of desired variables
674 ! found in Vtable.
675 ! ----
676 MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
677 ! maxvar is defined in table.mod
679 if ( gfld%discipline .eq. g2code(1,i) .and. !Discipline
680 & gfld%ipdtmpl(1) .eq. g2code(2,i) .and. !Category
681 & gfld%ipdtmpl(2) .eq. g2code(3,i) .and. !Parameter
682 & gfld%ipdtmpl(10) .eq. g2code(4,i) .and. !Elevation
683 & gfld%ipdtnum .eq. g2code(5,i)) then !Template
685 call gf_free(gfld)
686 call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
687 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
688 & gfld%ipdtmpl(2))
690 !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
691 my_field=namvar(i)
693 !MGD if (debug_level .gt. 50) then
694 !MGD write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
695 !MGD write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
696 !MGD endif
697 ! The following if-block is commented out since equivalent info can be obtained from g2print
698 ! if (debug_level .gt. 1000) then
699 ! fldmax=gfld%fld(1)
700 ! fldmin=gfld%fld(1)
701 ! sum=gfld%fld(1)
702 ! do j=2,gfld%ndpts
703 ! if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
704 ! if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
705 ! sum=sum+gfld%fld(j)
706 ! enddo ! gfld%ndpts
707 ! call mprintf(.true.,DEBUG,'G2 FIELD=%s MIN=%f AVG=%f MAX=%f',
708 ! & newline=.true., s1=pabbrev, f1=fldmin, f2=sum/gfld%ndpts,
709 ! & f3=fldmax)
710 ! endif
712 ! need to match up soil levels with those requested.
713 ! For the Vtable levels, -88 = all levels, -99 = missing. The units
714 ! vary depending on the level code (e.g. 106 = cm, 103 = m).
715 ! The grib2 standard allows scaling of the units, so make sure the soil level
716 ! units are in cm (as used in the Vtable).
717 if ( gfld%ipdtmpl(10) .eq. 106 ) then
718 if ( ( gfld%ipdtmpl(14) .EQ. -1*(2**07-1) ) .AND.
719 ! & ( gfld%ipdtmpl(15) .EQ. -1*(2**31-1) ) ) THEN ! Some compilers cannot
720 ! handle the initial 2**31
721 ! part of the computation,
722 ! which is an arithmetic
723 ! overflow on 32 bit signed ints
724 & ( gfld%ipdtmpl(15) .EQ. -2147483647 ) ) THEN
725 ! special UM grib2
726 glevel1 = gfld%ipdtmpl(12)
727 glevel2 = gfld%ipdtmpl(11)
728 else
729 glevel1 = 100. * gfld%ipdtmpl(12)*
730 & (10.**(-1.*gfld%ipdtmpl(11)))
731 glevel2 = 100. * gfld%ipdtmpl(15)*
732 & (10.**(-1.*gfld%ipdtmpl(14)))
733 end if
734 TMP8LOOP: do j = 1, maxvar
735 if ((g2code(4,j) .eq. 106) .and.
736 & (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
737 & (glevel1 .eq. level1(j)) .and.
738 & ((glevel2 .eq. level2(j)) .or.
739 & (level2(j) .le. -88))) then
740 my_field = namvar(j)
741 exit TMP8LOOP
742 endif
743 enddo TMP8LOOP
744 if (j .gt. maxvar ) then
745 write(6,'(a,i6,a)') 'Subsoil level ',
746 & gfld%ipdtmpl(12),
747 & ' in the GRIB2 file, was not found in the Vtable'
748 cycle MATCH_LOOP
749 endif
750 !MGD if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
751 endif
753 ! Level (eg. 10000 mb)
754 if(gfld%ipdtmpl(10).eq.100) then
755 ! Pressure level (range from 1000mb to 0mb)
756 level=gfld%ipdtmpl(12) *
757 & (10. ** (-1. * gfld%ipdtmpl(11)))
758 if ( level .lt. pmin ) cycle MATCH_LOOP
759 elseif((gfld%ipdtmpl(10).eq.105).or.
760 & (gfld%ipdtmpl(10).eq.118))then
761 ! Hybrid level (range from 1 to N)
762 level=gfld%ipdtmpl(12)
763 elseif(gfld%ipdtmpl(10).eq.104) then
764 ! Sigma level (range from 10000 to 0)
765 level=gfld%ipdtmpl(12)
766 elseif(gfld%ipdtmpl(10).eq.101) then
767 ! MSL
768 level=201300.
769 elseif(gfld%ipdtmpl(10).eq.103) then
770 ! Height above ground (m)
771 if (gfld%ipdtmpl(12) .eq. 2. .or.
772 & gfld%ipdtmpl(12) .eq. 1000. .or. ! temp fix for hrrr maxref
773 & gfld%ipdtmpl(12) .eq. 10. ) then
774 level=200100.
775 else
776 cycle MATCH_LOOP
777 endif
778 elseif((gfld%ipdtmpl(10).ge.206 .and.
779 & gfld%ipdtmpl(10).le.234) .or.
780 & (gfld%ipdtmpl(10).ge.242 .and.
781 & gfld%ipdtmpl(10).le.254) .or.
782 & (gfld%ipdtmpl(10).eq.200) .or.
783 & (gfld%ipdtmpl(10).eq.10) ) then
784 ! NCEP cloud layers used for plotting
785 level=200100.
786 elseif(gfld%ipdtmpl(10).eq.106.or.
787 & gfld%ipdtmpl(10).eq.1) then
788 ! Misc near ground/surface levels
789 level=200100.
790 elseif(gfld%ipdtmpl(10).eq.6) then
791 ! Level of Max wind
792 level=200100.
793 elseif(gfld%ipdtmpl(10).eq.7) then
794 ! Tropopause
795 level=200100.
796 else
797 ! If we are here then the Vtable contains a level code
798 ! which we cannot handle. Write an info message and skip it.
799 call mprintf(.true.,INFORM,"Rd_grib2 does not know abou
800 &t level code %i (field = %s). Skipping this field. If you want thi
801 &s level, rd_grib2.F must be modified", i1 = gfld%ipdtmpl(10),
802 & s1 = my_field )
803 cycle MATCH_LOOP
804 endif
805 iplvl = int(level)
807 ! Store the unpacked 2D slab from the GRIB2 record
808 allocate(hold_array(gfld%ngrdpts))
809 do j=1,gfld%ngrdpts
810 hold_array(j)=gfld%fld(j)
811 enddo
813 ! Some grids need to be reordered. Until we get an example, this is
814 ! a placeholder
815 ! call reorder_it (hold_array, map%nx, map%ny, map%dx,
816 ! & map%dy, iorder)
818 ! When we have reached this point, we have a data array ARRAY
819 ! which has some data we want to save, with field name FIELD
820 ! at pressure level LEVEL (Pa). Dimensions of this data are
821 ! map%nx and map%ny. Put this data into storage.
823 !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
824 !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
825 ! call mprintf(.true.,DEBUG,"Calling put_storage for
826 ! &level = %i , field = %s , g2level = %i ", newline=.true.,
827 ! & i1 = iplvl, s1 = my_field, i2 = gfld%ipdtmpl(12) )
829 call put_storage(iplvl,my_field,
830 & reshape(hold_array(1:map%nx*map%ny),
831 & (/map%nx, map%ny/)), map%nx,map%ny)
832 deallocate(hold_array)
834 ! If Specific Humidity is present on hybrid levels AND
835 ! upper-air RH is missing, see if we can compute RH from
836 ! Specific Humidity.
837 if (.not. is_there(iplvl, 'RH') .and.
838 & is_there(iplvl, 'SH') .and.
839 & is_there(iplvl, 'TT') .and.
840 & is_there(iplvl, 'P')) then
841 call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
842 !call llstor_remove(iplvl, 'SH') !We are done with SH
843 endif
845 ! If Specific Humidity is present on hybrid levels AND
846 ! upper-air RH is missing, see if we can compute RH from
847 ! Specific Humidity - v2
848 if (.not. is_there(iplvl, 'RH') .and.
849 & is_there(iplvl, 'SPECHUMD') .and.
850 & is_there(iplvl, 'THETA') .and.
851 & is_there(iplvl, 'TT')) then
852 call g2_compute_rh_spechumd_upa2(map%nx,map%ny,iplvl)
853 endif
855 ! If Temperature and Theta are present on hybrid levels AND
856 ! upper-air PRESSURE is missing, see if we can compute PRESSURE from
857 ! Temperature and Theta
858 if (.not. is_there(iplvl, 'PRESSURE') .and.
859 & is_there(iplvl, 'THETA') .and.
860 & is_there(iplvl, 'TT')) then
861 call g2_compute_pressure_tth_upa(map%nx,map%ny,iplvl)
862 endif
864 ith=ith+1
865 exit MATCH_LOOP
867 endif ! Selected param.
870 enddo MATCH_LOOP
872 ! Deallocate arrays decoding GRIB2 record.
873 call gf_free(gfld)
875 enddo NUM_FIELDS
878 enddo VERSION ! skgb
881 if ( debug_level .gt. 100 ) then
882 call mprintf (.true.,DEBUG,
883 & "G2 total number of fields found = %i ",newline=.true.,i1=itot)
884 end if
886 CALL BACLOSE(junit,IOS)
888 nullify(gfld%local) ! must be nullified before opening next file
889 ireaderr=1
890 else
891 call mprintf (.true.,DEBUG,"open status failed because %i ",
892 & newline=.true., i1=ios)
893 hdate = '9999-99-99_99:99:99'
894 ireaderr=2
895 endif ! ireaderr check
897 END subroutine rd_grib2
899 !*****************************************************************************!
900 ! Subroutine edition_num !
901 ! !
902 ! Purpose: !
903 ! Read one record from the input GRIB2 file. Based on the information in !
904 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
905 ! the GRIB2 record is one to process or to skip. If the field is one we !
906 ! want to keep, extract the data from the GRIB2 record, and pass the data !
907 ! back to the calling routine. !
908 ! !
909 ! Argument list: !
910 ! Input: !
911 ! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
912 ! unit number, since we do not do Fortran I/O for the GRIB2 !
913 ! files. Nor is it a UNIX File Descriptor returned from a C !
914 ! OPEN statement. It is really just an array index to the !
915 ! array (IUARR) where the UNIX File Descriptor values are !
916 ! stored. !
917 ! GRIB2FILE: File name to open, if it is not already open. !
918 ! !
919 ! Output: !
920 ! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
921 ! IERR : Error flag: 0 - no error on read from GRIB2 file. !
922 ! 1 - Hit the end of the GRIB2 file. !
923 ! 2 - The file GRIBFLNM we tried to open does !
924 ! not exist. !
925 ! Author: Paula McCaslin !
926 ! NOAA/FSL !
927 ! Sept 2004 !
928 !*****************************************************************************!
930 SUBROUTINE edition_num(junit, gribflnm,
931 & grib_edition, ireaderr)
933 use grib_mod
934 use params
935 use module_debug
937 parameter(msk1=32000,msk2=4000)
938 character(len=1),allocatable,dimension(:) :: cgrib
939 integer :: listsec0(3)
940 integer :: listsec1(13)
941 character(len=*) :: gribflnm
942 integer :: lskip, lgrib
943 integer :: junit
944 integer :: grib_edition
945 integer :: i, j, ireaderr
946 integer :: currlen
948 character(len=4) :: ctemp
949 character(len=4),parameter :: grib='GRIB',c7777='7777'
951 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
952 C SET ARGUMENTS
954 itot=0
955 icount=0
956 iseek=0
957 lskip=0
958 lgrib=0
959 currlen=0
961 !/* IOS Return Codes from BACIO: */
962 !/* 0 All was well */
963 !/* -1 Tried to open read only _and_ write only */
964 !/* -2 Tried to read and write in the same call */
965 !/* -3 Internal failure in name processing */
966 !/* -4 Failure in opening file */
967 !/* -5 Tried to read on a write-only file */
968 !/* -6 Failed in read to find the 'start' location */
969 !/* -7 Tried to write to a read only file */
970 !/* -8 Failed in write to find the 'start' location */
971 !/* -9 Error in close */
972 !/* -10 Read or wrote fewer data than requested */
974 !if ireaderr =1 we have hit the end of a file.
975 !if ireaderr =2 we have hit the end of all the files.
976 !if ireaderr =3 beginning characters 'GRIB' not found
978 ! Open a byte-addressable file.
979 CALL BAOPENR(junit,gribflnm,IOS)
980 if (ios.eq.0) then
982 ! Search opend file for the next GRIB2 messege (record).
983 call skgb(junit,iseek,msk1,lskip,lgrib)
985 ! Check for EOF, or problem
986 call mprintf((lgrib.eq.0),ERROR,
987 & "Grib2 file or date problem, stopping in edition_num.",
988 & newline=.true.)
990 ! Check size, if needed allocate more memory.
991 if (lgrib.gt.currlen) then
992 if (allocated(cgrib)) deallocate(cgrib)
993 allocate(cgrib(lgrib),stat=is)
994 currlen=lgrib
995 endif
997 ! Read a given number of bytes from unblocked file.
998 call baread(junit,lskip,lgrib,lengrib,cgrib)
1000 ! Check for beginning of GRIB message in the first 100 bytes
1001 istart=0
1002 do j=1,100
1003 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
1004 if (ctemp.eq.grib ) then
1005 istart=j
1006 exit
1007 endif
1008 enddo
1009 if (istart.eq.0) then
1010 ireaderr=3
1011 print*, "The beginning 4 characters >GRIB< were not found."
1012 endif
1014 ! Unpack Section 0 - Indicator Section to extract GRIB edition field
1015 iofst=8*(istart+5)
1016 call gbyte(cgrib,discipline,iofst,8) ! Discipline
1017 iofst=iofst+8
1018 call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
1020 ! print *, 'ungrib - grib edition num', grib_edition
1021 CALL BACLOSE(junit,IOS)
1022 ireaderr=1
1023 else if (ios .eq. -4) then
1024 call mprintf(.true.,ERROR,
1025 & "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
1026 else
1027 print *,'edition_num: open status failed because',ios,gribflnm
1028 ireaderr=2
1029 endif ! ireaderr check
1031 END subroutine edition_num
1033 !*****************************************************************************!
1035 SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
1036 ! Compute relative humidity from specific humidity in the upper air.
1037 use storage_module
1038 implicit none
1039 integer :: ix, jx
1040 integer :: iiplvl
1041 real :: lat1, lon1, dx, dy
1042 real, dimension(ix,jx) :: T, P, RH, Q
1044 real, parameter :: svp1=611.2
1045 real, parameter :: svp2=17.67
1046 real, parameter :: svp3=29.65
1047 real, parameter :: svpt0=273.15
1048 real, parameter :: eps = 0.622
1050 real startlat, startlon, deltalat, deltalon
1052 call get_storage(iiplvl, 'P', P, ix, jx)
1053 call get_storage(iiplvl, 'TT', T, ix, jx)
1054 call get_storage(iiplvl, 'SH', Q, ix, jx)
1056 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1058 call put_storage(iiplvl, 'RH', rh, ix, jx)
1060 end subroutine g2_compute_rh_spechumd_upa
1062 !*****************************************************************************!
1064 SUBROUTINE g2_compute_rh_spechumd_upa2(ix, jx, iiplvl)
1065 ! Compute relative humidity from specific humidity in the upper air.
1066 use storage_module
1067 implicit none
1068 integer :: ix, jx
1069 integer :: iiplvl
1070 real :: lat1, lon1, dx, dy
1071 real, dimension(ix,jx) :: T, TH, RH, Q, P
1073 real, parameter :: svp1=611.2
1074 real, parameter :: svp2=17.67
1075 real, parameter :: svp3=29.65
1076 real, parameter :: svpt0=273.15
1077 real, parameter :: eps = 0.622
1079 real startlat, startlon, deltalat, deltalon
1081 call get_storage(iiplvl, 'THETA', TH, ix, jx)
1082 call get_storage(iiplvl, 'TT', T, ix, jx)
1083 call get_storage(iiplvl, 'SPECHUMD', Q, ix, jx)
1085 p=1.e5*(t/th)**(1005/287.05)
1087 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1089 call put_storage(iiplvl, 'RH', rh, ix, jx)
1091 end subroutine g2_compute_rh_spechumd_upa2
1093 !*****************************************************************************!
1095 SUBROUTINE g2_compute_pressure_tth_upa(ix, jx, iiplvl)
1096 ! Compute relative humidity from specific humidity in the upper air.
1097 use storage_module
1098 implicit none
1099 integer :: ix, jx
1100 integer :: iiplvl
1101 real :: lat1, lon1, dx, dy
1102 real, dimension(ix,jx) :: T, TH, P
1104 real, parameter :: svp1=611.2
1105 real, parameter :: svp2=17.67
1106 real, parameter :: svp3=29.65
1107 real, parameter :: svpt0=273.15
1108 real, parameter :: eps = 0.622
1110 real startlat, startlon, deltalat, deltalon
1112 call get_storage(iiplvl, 'THETA', TH, ix, jx)
1113 call get_storage(iiplvl, 'TT', T, ix, jx)
1115 p=1.e5*(t/th)**(1005/287.05)
1117 call put_storage(iiplvl, 'PRESSURE', p, ix, jx)
1119 end subroutine g2_compute_pressure_tth_upa
1121 !*****************************************************************************!
1123 subroutine ncep_grid_num (pnum)
1124 !
1125 ! Find the grib number for descriptive labelling.
1126 ! Grib2 doesn't have a grid-number entry, so we have to figure it out
1127 ! from the parameters.
1128 !
1129 use gridinfo ! Included to define map%
1130 integer :: pnum
1131 real, parameter :: eps = .01
1132 character (len=8) :: tmp8
1134 ! write(6,*) 'begin ncep_grid_num'
1135 ! write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
1136 tmp8 = ' '
1137 if (pnum .eq. 30) then ! lambert conformal
1138 if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
1139 write(tmp8,'("GRID 218")')
1140 else if (abs(map%dx - 40.63525) .lt. eps
1141 & .and. map%nx .eq. 185) then
1142 write(tmp8,'("GRID 212")')
1143 else if (abs(map%dx - 40.63525) .lt. eps
1144 & .and. map%nx .eq. 151) then
1145 write(tmp8,'("GRID 236")')
1146 else if (abs(map%dx - 81.2705) .lt. eps
1147 & .and. map%nx .eq. 93) then
1148 write(tmp8,'("GRID 211")')
1149 else if (abs (map%dx - 32.46341) .lt. eps
1150 & .and. map%nx .eq. 349) then
1151 write(tmp8,'("GRID 221")')
1152 else if (abs(map%dx - 20.317625) .lt. eps
1153 & .and. map%nx .eq. 301) then
1154 write(tmp8,'("GRID 252")')
1155 else if (abs(map%dx - 13.545087) .lt. eps
1156 & .and. map%nx .eq. 451) then
1157 write(tmp8,'("GRID 130")')
1158 endif
1159 else if (pnum .eq. 20) then ! polar stereographic
1160 if (abs(map%dx - 15.0) .lt. eps) then
1161 write(tmp8,'("GRID 88")')
1162 endif
1163 else if (pnum .eq. 0) then ! lat/lon
1164 if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
1165 write(tmp8,'("GRID 3")')
1166 else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
1167 write(tmp8,'("GRID 4")')
1168 endif
1169 endif
1170 map%source(25:32) = tmp8
1171 ! write(6,*) 'map%source = ',map%source
1172 end subroutine ncep_grid_num
1173 !*****************************************************************************!
1175 function earth_radius (icode, iscale, irad_m)
1176 ! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
1177 use module_debug
1178 real :: earth_radius
1179 integer :: icode
1180 integer :: iscale, irad_m
1181 if ( icode .eq. 0 ) then
1182 earth_radius = 6367470. * .001
1183 else if ( icode .eq. 1) then
1184 earth_radius = 0.001 * float(irad_m) / 10**iscale
1185 else if ( icode .eq. 6 ) then
1186 earth_radius = 6371229. * .001
1187 else if ( icode .eq. 8 ) then
1188 earth_radius = 6371200. * .001
1189 else
1190 call mprintf(.true.,ERROR,
1191 & "unknown earth radius for code %i",newline=.true.,i1=icode)
1192 endif
1193 end function earth_radius