| blob | 5a584384c7049f5189c48fcc0eb2dbb6a7a9121c |
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 if (icenter .eq. 78 .or. icenter .eq. 79 ) then
330 map%source = 'DWD'
331 else
332 map%source = 'unknown model and orig center'
333 end if
334 if (icenter.eq.7) then
335 if (iprocess.eq.81 .or. iprocess.eq.82 .or.
336 & iprocess.eq.96) then
337 ! pmin should not be 15 or 40 hPa for gfs files. Use the next highest level.
338 if (pmin .eq. 1500.) then
339 pmin = 1000.
340 else if (pmin .eq. 4000.) then
341 pmin = 3000.
342 endif
343 endif
344 endif
346 call mprintf(.true.,DEBUG,"G2 source is = %s ",
347 & newline=.true., s1=map%source)
349 !--
351 ! Store information about the grid containing the data.
352 ! This stuff gets stored in the MAP variable, as defined in
353 ! module GRIDINFO.
355 map%startloc = 'SWCORNER'
356 map%grid_wind = .true.
358 if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
359 map%igrid = 0
360 map%nx = gfld%igdtmpl(8)
361 map%ny = gfld%igdtmpl(9)
362 map%dx = gfld%igdtmpl(17)
363 map%dy = gfld%igdtmpl(18)
364 map%lat1 = gfld%igdtmpl(12)
365 map%lon1 = gfld%igdtmpl(13)
366 write(tmp8,'(b8.8)') gfld%igdtmpl(14)
367 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
368 map%r_earth = earth_radius (gfld%igdtmpl(1),
369 & gfld%igdtmpl(2),gfld%igdtmpl(3))
371 ! Fix for NCEP 1/12 degree grids (e.g. rtgsst)
372 if (icenter .eq. 7 .and. map%dx .eq. 83000. .and. map%nx
373 & .eq. 4320) then
374 map%lat1 = 89958333.
375 map%lon1 = 41667.
376 map%dx = 83333.333 * sign(1.0,map%dx)
377 map%dy = 83333.333 * sign(1.0,map%dy)
378 endif
380 if ((gfld%igdtmpl(10) .eq. 0).OR.
381 & (gfld%igdtmpl(10) .eq. 255)) THEN
382 ! Scale lat/lon values to 0-180, default range is 1e6.
383 map%lat1 = map%lat1/scale_factor
384 map%lon1 = map%lon1/scale_factor
385 ! Scale dx/dy values to degrees, default range is 1e6.
386 map%dx = map%dx/scale_factor
387 map%dy = map%dy/scale_factor
388 else
389 ! Basic angle and subdivisions are non-zero (not tested)
390 map%lat1 = map%lat1 *
391 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
392 map%lon1 = map%lon1 *
393 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
394 map%dx = map%dx *
395 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
396 map%dy = map%dy *
397 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
398 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
399 &has not been tested, continuing anyway")
400 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
401 & has not been tested, continuing anyway")
402 endif
405 ! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
406 ! In WPS, the sign of dy indicates the direction of the scan.
407 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
408 read(tmp8,'(1x,i1)') jscan
409 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
410 map%dy = -1. * map%dy
411 endif
412 ! if ( map%lat1 .gt. gfld%igdtmpl(15) .and.
413 ! & map%dy .gt. 0. ) then
414 ! map%dy = -1. * map%dy
415 ! write(6,*) 'Resetting map%dy for iprocess = ',iprocess
416 ! endif
418 elseif (gfld%igdtnum.eq.10) then ! Mercator Grid.
419 map%igrid = 1
420 map%nx = gfld%igdtmpl(8)
421 map%ny = gfld%igdtmpl(9)
422 map%lov = 0.
423 map%truelat1 = gfld%igdtmpl(13) / scale_factor
424 map%truelat2 = 0.
425 map%dx = gfld%igdtmpl(18) / scale_factor
426 map%dy = gfld%igdtmpl(19) / scale_factor
427 map%lat1 = gfld%igdtmpl(10) / scale_factor
428 map%lon1 = gfld%igdtmpl(11) / scale_factor
429 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
430 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
431 map%r_earth = earth_radius (gfld%igdtmpl(1),
432 & gfld%igdtmpl(2),gfld%igdtmpl(3))
434 elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
435 map%igrid = 5
436 map%nx = gfld%igdtmpl(8)
437 map%ny = gfld%igdtmpl(9)
438 map%lov = gfld%igdtmpl(14) / scale_factor
439 map%truelat1 = 60.
440 map%truelat2 = 91.
441 map%dx = gfld%igdtmpl(15) / scale_factor
442 map%dy = gfld%igdtmpl(16) / scale_factor
443 map%lat1 = gfld%igdtmpl(10) / scale_factor
444 map%lon1 = gfld%igdtmpl(11) / scale_factor
445 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
446 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
447 map%r_earth = earth_radius (gfld%igdtmpl(1),
448 & gfld%igdtmpl(2),gfld%igdtmpl(3))
450 elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
451 map%igrid = 3
452 map%nx = gfld%igdtmpl(8)
453 map%ny = gfld%igdtmpl(9)
454 map%lov = gfld%igdtmpl(14) / scale_factor
455 map%truelat1 = gfld%igdtmpl(19) / scale_factor
456 map%truelat2 = gfld%igdtmpl(20) / scale_factor
457 map%dx = gfld%igdtmpl(15) / scale_factor
458 map%dy = gfld%igdtmpl(16) / scale_factor
459 map%lat1 = gfld%igdtmpl(10) / scale_factor
460 map%lon1 = gfld%igdtmpl(11) / scale_factor
461 write(tmp8,'(b8.8)') gfld%igdtmpl(12)
462 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
463 map%r_earth = earth_radius (gfld%igdtmpl(1),
464 & gfld%igdtmpl(2),gfld%igdtmpl(3))
466 elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
467 map%igrid = 4
468 map%nx = gfld%igdtmpl(8) ! Ni - # of points along a parallel
469 map%ny = gfld%igdtmpl(9) ! Nj - # of points along meridian
470 map%dx = gfld%igdtmpl(17) ! Di - i direction increment
471 map%dy = gfld%igdtmpl(18) ! N - # of parallels between pole and equator
472 map%lat1 = gfld%igdtmpl(12) ! La1 - lat of 1st grid point
473 map%lon1 = gfld%igdtmpl(13) ! Lo1 - lon of 1st grid point
474 write(tmp8,'(b8.8)') gfld%igdtmpl(14) ! resolution/component flag
475 if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
476 map%r_earth = earth_radius (gfld%igdtmpl(1),
477 & gfld%igdtmpl(2),gfld%igdtmpl(3))
479 ! Scale dx/dy values to degrees, default range is 1e6.
480 if (map%dx.gt.10000) then
481 map%dx = map%dx/scale_factor
482 endif
483 if (map%dy.gt.10000) then
484 map%dy = (map%dy/scale_factor)*(-1)
485 endif
487 ! Fix for zonal shift in CFSR data, following a similar fix
488 ! for global lat-lon data in rd_grib1.F
489 if ( ABS(map%nx * map%dx - 360.0) < 1.0 ) then
490 if (ABS(map%dx - (360./real(map%nx))) > 0.00001) then
491 write(0,*) 'CFSR fix: recomputing delta-longitude'
492 map%dx = 360./real(map%nx)
493 endif
494 endif
496 ! Scale lat/lon values to 0-180, default range is 1e6.
497 if (abs(map%lat1).ge.scale_factor) then
498 map%lat1 = map%lat1/scale_factor
499 endif
500 if (map%lon1.ge.scale_factor) then
501 map%lon1 = map%lon1/scale_factor
502 endif
503 if ( debug_level .gt. 2 ) then
504 call mprintf(.true.,DEBUG,
505 & "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
506 & newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
507 & i1=nint(map%dy))
508 end if
509 elseif (gfld%igdtnum.eq.32769) then ! Arakawa Non-E Staggered grid.
510 map%igrid = 6
511 map%nx = gfld%igdtmpl(8)
512 map%ny = gfld%igdtmpl(9)
513 map%dx = gfld%igdtmpl(17)
514 map%dy = gfld%igdtmpl(18)
515 map%lat1 = gfld%igdtmpl(12)
516 map%lon1 = gfld%igdtmpl(13)
517 map%lat0 = gfld%igdtmpl(15)
518 map%lon0 = gfld%igdtmpl(16)
519 map%r_earth = earth_radius (gfld%igdtmpl(1),
520 & gfld%igdtmpl(2),gfld%igdtmpl(3))
521 if ((gfld%igdtmpl(10) .eq. 0).OR.
522 & (gfld%igdtmpl(10) .eq. 255)) THEN
523 map%lat1 = map%lat1/scale_factor
524 map%lon1 = map%lon1/scale_factor
525 map%lat0 = map%lat0/scale_factor
526 map%lon0 = map%lon0/scale_factor
527 map%dx = map%dx/scale_factor/1.e3
528 map%dy = map%dy/scale_factor/1.e3
529 else
530 ! Basic angle and subdivisions are non-zero (not tested)
531 map%lat1 = map%lat1 *
532 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
533 map%lon1 = map%lon1 *
534 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
535 map%dx = map%dx *
536 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
537 map%dy = map%dy *
538 & (gfld%igdtmpl(11)/gfld%igdtmpl(10))
539 call mprintf(.true.,STDOUT,"WARNING - Basic angle option
540 &has not been tested, continuing anyway")
541 call mprintf(.true.,LOGFILE,"WARNING - Basic angle option
542 & has not been tested, continuing anyway")
543 endif
545 else
546 call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
547 & newline=.true.,i1=gfld%igdtnum)
548 call mprintf(.true.,STDOUT,
549 & "ungrib understands projections 0, 20, 30, and 40",
550 & newline=.true.)
551 call mprintf(.true.,LOGFILE,
552 & "GRIB2 Unknown Projection: %i",
553 & newline=.true.,i1=gfld%igdtnum)
554 call mprintf(.true.,LOGFILE,
555 & "ungrib understands projections 0, 10, 20, 30, and 40",
556 & newline=.true.)
557 ! If the projection is not known, then it can't be processed by metgrid/plotfmt
558 stop 'Stop in rd_grib2'
559 endif
561 call mprintf(.true.,DEBUG,"G2 igrid = %i , dx = %f , dy = %
562 &f ", newline=.true., i1 = map%igrid, f1=map%dx, f2=map%dy)
564 if (icenter.eq.7) then
565 call ncep_grid_num (gfld%igdtnum)
566 endif
568 ! Deallocate arrays decoding GRIB2 record.
569 call gf_free(gfld)
571 endif ! "first" if-block
573 ! ----
575 ! Continue to unpack GRIB2 field.
576 NUM_FIELDS: do n = 1, numfields
577 ! e.g. U and V would =2, otherwise its usually =1
578 call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
579 if (ierr.ne.0) then
580 write(*,*) ' ERROR extracting field gf_getfld = ',ierr
581 cycle
582 endif
584 ! The JMA GSM has two different grids in the same GRIB file, so we need
585 ! to process the map info for each field separately. If any other centers do
586 ! this, then processing will need to be added here, too.
588 if (icenter .eq. 34 .and. gfld%igdtnum.eq.0) then
589 map%nx = gfld%igdtmpl(8)
590 map%ny = gfld%igdtmpl(9)
591 map%dx = gfld%igdtmpl(17)
592 map%dy = gfld%igdtmpl(18)
593 ! Scale dx/dy values to degrees, default range is 1e6.
594 if (map%dx.gt.10000) then
595 map%dx = map%dx/scale_factor
596 endif
597 if (map%dy.gt.10000) then
598 map%dy = map%dy/scale_factor
599 endif
600 write(tmp8,'(b8.8)') gfld%igdtmpl(19)
601 read(tmp8,'(1x,i1)') jscan
602 write(0,*) 'gfld%igdtmpl(19) = ',gfld%igdtmpl(19),
603 & ' jscan = ',jscan
604 if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
605 map%dy = -1. * map%dy
606 endif
607 endif ! JMA spectral
609 ! ------------------------------------
610 ! Additional print information for developer.
611 if ( debug_level .GT. 1000 ) then
612 !MGD print *
613 !MGD print *,'G2 FIELD ',n
614 !MGD if (n==1) then
615 !MGD print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
616 !MGD print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
617 !MGD endif
618 !MGD if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
619 !MGD print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
620 !MGD endif
621 !MGD print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
622 !MGD & gfld%numoct_opt,gfld%interp_opt,
623 !MGD & gfld%igdtnum
624 !MGD print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
625 !MGD & (gfld%igdtmpl(j),j=1,gfld%igdtlen)
626 !MGD if ( gfld%num_opt .eq. 0 ) then
627 !MGD print *,'G2 NO Section 3 List Defining No. of Data Points.'
628 !MGD else
629 !MGD print *,'G2 Section 3 Optional List: ',
630 !MGD & (gfld%list_opt(j),j=1,gfld%num_opt)
631 !MGD endif
632 !MGD print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
633 !MGD & (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
635 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
636 & gfld%ipdtmpl(2))
637 !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
638 !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
639 !MGD print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
641 !MGD if ( gfld%num_coord .eq. 0 ) then
642 !MGD print *,'G2 NO Optional Vertical Coordinate List.'
643 !MGD else
644 !MGD print *,'G2 Section 4 Optional Coordinates: ',
645 !MGD & (gfld%coord_list(j),j=1,gfld%num_coord)
646 !MGD endif
647 if ( gfld%ibmap .ne. 255 ) then
648 call mprintf(.true.,DEBUG,
649 & 'G2 Num. of Data Points = %i with BIT-MAP %i',
650 & newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
651 else
652 call mprintf(.true.,DEBUG,
653 & 'G2 Num. of Data Points = %i NO BIT-MAP',
654 & newline=.true., i1=gfld%ndpts)
655 endif
656 !MGD print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
657 !MGD & (gfld%idrtmpl(j),j=1,gfld%idrtlen)
658 endif ! Additional Print information
659 ! ------------------------------------
661 ! do i = 1, maxvar
662 ! write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
663 ! enddo
664 !MGD if (debug_level .gt. 50) then
665 !MGD write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
666 !MGD & gfld%ipdtmpl(2),gfld%ipdtmpl(10)
667 !MGD endif
668 call mprintf(.true.,DEBUG,"G2 Searching the g2code array (Vta
669 &ble) for this grib field %i %i %i %i %i %i ", newline=.true.,
670 & i1 = gfld%discipline, i2 = gfld%ipdtmpl(1),
671 & i3 = gfld%ipdtmpl(2), i4 = gfld%ipdtmpl(10),
672 & i5 = gfld%ipdtmpl(12), i6 = gfld%ipdtnum )
675 ! Test this data record against list of desired variables
676 ! found in Vtable.
677 ! ----
678 MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
679 ! maxvar is defined in table.mod
681 if ( gfld%discipline .eq. g2code(1,i) .and. !Discipline
682 & gfld%ipdtmpl(1) .eq. g2code(2,i) .and. !Category
683 & gfld%ipdtmpl(2) .eq. g2code(3,i) .and. !Parameter
684 & gfld%ipdtmpl(10) .eq. g2code(4,i) .and. !Elevation
685 & gfld%ipdtnum .eq. g2code(5,i)) then !Template
687 call gf_free(gfld)
688 call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
689 pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
690 & gfld%ipdtmpl(2))
692 !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
693 my_field=namvar(i)
695 !MGD if (debug_level .gt. 50) then
696 !MGD write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
697 !MGD write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
698 !MGD endif
699 ! The following if-block is commented out since equivalent info can be obtained from g2print
700 ! if (debug_level .gt. 1000) then
701 ! fldmax=gfld%fld(1)
702 ! fldmin=gfld%fld(1)
703 ! sum=gfld%fld(1)
704 ! do j=2,gfld%ndpts
705 ! if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
706 ! if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
707 ! sum=sum+gfld%fld(j)
708 ! enddo ! gfld%ndpts
709 ! call mprintf(.true.,DEBUG,'G2 FIELD=%s MIN=%f AVG=%f MAX=%f',
710 ! & newline=.true., s1=pabbrev, f1=fldmin, f2=sum/gfld%ndpts,
711 ! & f3=fldmax)
712 ! endif
714 ! need to match up soil levels with those requested.
715 ! For the Vtable levels, -88 = all levels, -99 = missing. The units
716 ! vary depending on the level code (e.g. 106 = cm, 103 = m).
717 ! The grib2 standard allows scaling of the units, so make sure the soil level
718 ! units are in cm (as used in the Vtable).
719 if ( gfld%ipdtmpl(10) .eq. 106 ) then
720 if ( ( gfld%ipdtmpl(14) .EQ. -1*(2**07-1) ) .AND.
721 ! & ( gfld%ipdtmpl(15) .EQ. -1*(2**31-1) ) ) THEN ! Some compilers cannot
722 ! handle the initial 2**31
723 ! part of the computation,
724 ! which is an arithmetic
725 ! overflow on 32 bit signed ints
726 & ( gfld%ipdtmpl(15) .EQ. -2147483647 ) ) THEN
727 ! special UM grib2
728 glevel1 = gfld%ipdtmpl(12)
729 glevel2 = gfld%ipdtmpl(11)
730 else
731 glevel1 = 100. * gfld%ipdtmpl(12)*
732 & (10.**(-1.*gfld%ipdtmpl(11)))
733 glevel2 = 100. * gfld%ipdtmpl(15)*
734 & (10.**(-1.*gfld%ipdtmpl(14)))
735 end if
736 TMP8LOOP: do j = 1, maxvar
737 if ((g2code(4,j) .eq. 106) .and.
738 & (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
739 & (glevel1 .eq. level1(j)) .and.
740 & ((glevel2 .eq. level2(j)) .or.
741 & (level2(j) .le. -88))) then
742 my_field = namvar(j)
743 exit TMP8LOOP
744 endif
745 enddo TMP8LOOP
746 if (j .gt. maxvar ) then
747 write(6,'(a,i6,a)') 'Subsoil level ',
748 & gfld%ipdtmpl(12),
749 & ' in the GRIB2 file, was not found in the Vtable'
750 cycle MATCH_LOOP
751 endif
752 !MGD if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
753 endif
755 ! Level (eg. 10000 mb)
756 if(gfld%ipdtmpl(10).eq.100) then
757 ! Pressure level (range from 1000mb to 0mb)
758 level=gfld%ipdtmpl(12) *
759 & (10. ** (-1. * gfld%ipdtmpl(11)))
760 if ( level .lt. pmin ) cycle MATCH_LOOP
761 elseif((gfld%ipdtmpl(10).eq.105).or.
762 & (gfld%ipdtmpl(10).eq.118).or.
763 & (gfld%ipdtmpl(10).eq.150))then
764 ! Hybrid level (range from 1 to N)
765 level=gfld%ipdtmpl(12)
766 elseif(gfld%ipdtmpl(10).eq.104) then
767 ! Sigma level (range from 10000 to 0)
768 level=gfld%ipdtmpl(12)
769 elseif(gfld%ipdtmpl(10).eq.101) then
770 ! MSL
771 level=201300.
772 elseif(gfld%ipdtmpl(10).eq.103) then
773 ! Height above ground (m)
774 if (gfld%ipdtmpl(12) .eq. 2. .or.
775 & gfld%ipdtmpl(12) .eq. 1000. .or. ! temp fix for hrrr maxref
776 & gfld%ipdtmpl(12) .eq. 10. ) then
777 level=200100.
778 else
779 cycle MATCH_LOOP
780 endif
781 elseif((gfld%ipdtmpl(10).ge.206 .and.
782 & gfld%ipdtmpl(10).le.234) .or.
783 & (gfld%ipdtmpl(10).ge.242 .and.
784 & gfld%ipdtmpl(10).le.254) .or.
785 & (gfld%ipdtmpl(10).eq.200) .or.
786 & (gfld%ipdtmpl(10).eq.10) ) then
787 ! NCEP cloud layers used for plotting
788 level=200100.
789 elseif(gfld%ipdtmpl(10).eq.106.or.
790 & gfld%ipdtmpl(10).eq.1) then
791 ! Misc near ground/surface levels
792 level=200100.
793 elseif(gfld%ipdtmpl(10).eq.6) then
794 ! Level of Max wind
795 level=200100.
796 elseif(gfld%ipdtmpl(10).eq.7) then
797 ! Tropopause
798 level=200100.
799 else
800 ! If we are here then the Vtable contains a level code
801 ! which we cannot handle. Write an info message and skip it.
802 call mprintf(.true.,INFORM,"Rd_grib2 does not know abou
803 &t level code %i (field = %s). Skipping this field. If you want thi
804 &s level, rd_grib2.F must be modified", i1 = gfld%ipdtmpl(10),
805 & s1 = my_field )
806 cycle MATCH_LOOP
807 endif
808 iplvl = int(level)
810 ! Store the unpacked 2D slab from the GRIB2 record
811 allocate(hold_array(gfld%ngrdpts))
812 do j=1,gfld%ngrdpts
813 hold_array(j)=gfld%fld(j)
814 enddo
816 ! Some grids need to be reordered. Until we get an example, this is
817 ! a placeholder
818 ! call reorder_it (hold_array, map%nx, map%ny, map%dx,
819 ! & map%dy, iorder)
821 ! When we have reached this point, we have a data array ARRAY
822 ! which has some data we want to save, with field name FIELD
823 ! at pressure level LEVEL (Pa). Dimensions of this data are
824 ! map%nx and map%ny. Put this data into storage.
826 !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
827 !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
828 ! call mprintf(.true.,DEBUG,"Calling put_storage for
829 ! &level = %i , field = %s , g2level = %i ", newline=.true.,
830 ! & i1 = iplvl, s1 = my_field, i2 = gfld%ipdtmpl(12) )
832 call put_storage(iplvl,my_field,
833 & reshape(hold_array(1:map%nx*map%ny),
834 & (/map%nx, map%ny/)), map%nx,map%ny)
835 deallocate(hold_array)
837 ! If Specific Humidity is present on hybrid levels AND
838 ! upper-air RH is missing, see if we can compute RH from
839 ! Specific Humidity.
840 if (.not. is_there(iplvl, 'RH') .and.
841 & is_there(iplvl, 'SH') .and.
842 & is_there(iplvl, 'TT') .and.
843 & is_there(iplvl, 'P')) then
844 call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
845 !call llstor_remove(iplvl, 'SH') !We are done with SH
846 endif
848 ! If Specific Humidity is present on hybrid levels AND
849 ! upper-air RH is missing, see if we can compute RH from
850 ! Specific Humidity - v2
851 if (.not. is_there(iplvl, 'RH') .and.
852 & is_there(iplvl, 'SPECHUMD') .and.
853 & is_there(iplvl, 'THETA') .and.
854 & is_there(iplvl, 'TT')) then
855 call g2_compute_rh_spechumd_upa2(map%nx,map%ny,iplvl)
856 endif
858 ! If Temperature and Theta are present on hybrid levels AND
859 ! upper-air PRESSURE is missing, see if we can compute PRESSURE from
860 ! Temperature and Theta
861 if (.not. is_there(iplvl, 'PRESSURE') .and.
862 & is_there(iplvl, 'THETA') .and.
863 & is_there(iplvl, 'TT')) then
864 call g2_compute_pressure_tth_upa(map%nx,map%ny,iplvl)
865 endif
867 ith=ith+1
868 exit MATCH_LOOP
870 endif ! Selected param.
873 enddo MATCH_LOOP
875 ! Deallocate arrays decoding GRIB2 record.
876 call gf_free(gfld)
878 enddo NUM_FIELDS
881 enddo VERSION ! skgb
884 if ( debug_level .gt. 100 ) then
885 call mprintf (.true.,DEBUG,
886 & "G2 total number of fields found = %i ",newline=.true.,i1=itot)
887 end if
889 CALL BACLOSE(junit,IOS)
891 nullify(gfld%local) ! must be nullified before opening next file
892 ireaderr=1
893 else
894 call mprintf (.true.,DEBUG,"open status failed because %i ",
895 & newline=.true., i1=ios)
896 hdate = '9999-99-99_99:99:99'
897 ireaderr=2
898 endif ! ireaderr check
900 END subroutine rd_grib2
902 !*****************************************************************************!
903 ! Subroutine edition_num !
904 ! !
905 ! Purpose: !
906 ! Read one record from the input GRIB2 file. Based on the information in !
907 ! the GRIB2 header and the user-defined Vtable, decide whether the field in!
908 ! the GRIB2 record is one to process or to skip. If the field is one we !
909 ! want to keep, extract the data from the GRIB2 record, and pass the data !
910 ! back to the calling routine. !
911 ! !
912 ! Argument list: !
913 ! Input: !
914 ! JUNIT : "Unit Number" to open and read from. Not really a Fortran !
915 ! unit number, since we do not do Fortran I/O for the GRIB2 !
916 ! files. Nor is it a UNIX File Descriptor returned from a C !
917 ! OPEN statement. It is really just an array index to the !
918 ! array (IUARR) where the UNIX File Descriptor values are !
919 ! stored. !
920 ! GRIB2FILE: File name to open, if it is not already open. !
921 ! !
922 ! Output: !
923 ! GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2 !
924 ! IERR : Error flag: 0 - no error on read from GRIB2 file. !
925 ! 1 - Hit the end of the GRIB2 file. !
926 ! 2 - The file GRIBFLNM we tried to open does !
927 ! not exist. !
928 ! Author: Paula McCaslin !
929 ! NOAA/FSL !
930 ! Sept 2004 !
931 !*****************************************************************************!
933 SUBROUTINE edition_num(junit, gribflnm,
934 & grib_edition, ireaderr)
936 use grib_mod
937 use params
938 use module_debug
940 parameter(msk1=32000,msk2=4000)
941 character(len=1),allocatable,dimension(:) :: cgrib
942 integer :: listsec0(3)
943 integer :: listsec1(13)
944 character(len=*) :: gribflnm
945 integer :: lskip, lgrib
946 integer :: junit
947 integer :: grib_edition
948 integer :: i, j, ireaderr
949 integer :: currlen
951 character(len=4) :: ctemp
952 character(len=4),parameter :: grib='GRIB',c7777='7777'
954 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
955 C SET ARGUMENTS
957 itot=0
958 icount=0
959 iseek=0
960 lskip=0
961 lgrib=0
962 currlen=0
964 !/* IOS Return Codes from BACIO: */
965 !/* 0 All was well */
966 !/* -1 Tried to open read only _and_ write only */
967 !/* -2 Tried to read and write in the same call */
968 !/* -3 Internal failure in name processing */
969 !/* -4 Failure in opening file */
970 !/* -5 Tried to read on a write-only file */
971 !/* -6 Failed in read to find the 'start' location */
972 !/* -7 Tried to write to a read only file */
973 !/* -8 Failed in write to find the 'start' location */
974 !/* -9 Error in close */
975 !/* -10 Read or wrote fewer data than requested */
977 !if ireaderr =1 we have hit the end of a file.
978 !if ireaderr =2 we have hit the end of all the files.
979 !if ireaderr =3 beginning characters 'GRIB' not found
981 ! Open a byte-addressable file.
982 CALL BAOPENR(junit,gribflnm,IOS)
983 if (ios.eq.0) then
985 ! Search opend file for the next GRIB2 messege (record).
986 call skgb(junit,iseek,msk1,lskip,lgrib)
988 ! Check for EOF, or problem
989 call mprintf((lgrib.eq.0),ERROR,
990 & "Grib2 file or date problem, stopping in edition_num.",
991 & newline=.true.)
993 ! Check size, if needed allocate more memory.
994 if (lgrib.gt.currlen) then
995 if (allocated(cgrib)) deallocate(cgrib)
996 allocate(cgrib(lgrib),stat=is)
997 currlen=lgrib
998 endif
1000 ! Read a given number of bytes from unblocked file.
1001 call baread(junit,lskip,lgrib,lengrib,cgrib)
1003 ! Check for beginning of GRIB message in the first 100 bytes
1004 istart=0
1005 do j=1,100
1006 ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
1007 if (ctemp.eq.grib ) then
1008 istart=j
1009 exit
1010 endif
1011 enddo
1012 if (istart.eq.0) then
1013 ireaderr=3
1014 print*, "The beginning 4 characters >GRIB< were not found."
1015 endif
1017 ! Unpack Section 0 - Indicator Section to extract GRIB edition field
1018 iofst=8*(istart+5)
1019 call gbyte(cgrib,discipline,iofst,8) ! Discipline
1020 iofst=iofst+8
1021 call gbyte(cgrib,grib_edition,iofst,8) ! GRIB edition number
1023 ! print *, 'ungrib - grib edition num', grib_edition
1024 CALL BACLOSE(junit,IOS)
1025 ireaderr=1
1026 else if (ios .eq. -4) then
1027 call mprintf(.true.,ERROR,
1028 & "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
1029 else
1030 print *,'edition_num: open status failed because',ios,gribflnm
1031 ireaderr=2
1032 endif ! ireaderr check
1034 END subroutine edition_num
1036 !*****************************************************************************!
1038 SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
1039 ! Compute relative humidity from specific humidity in the upper air.
1040 use storage_module
1041 implicit none
1042 integer :: ix, jx
1043 integer :: iiplvl
1044 real :: lat1, lon1, dx, dy
1045 real, dimension(ix,jx) :: T, P, RH, Q
1047 real, parameter :: svp1=611.2
1048 real, parameter :: svp2=17.67
1049 real, parameter :: svp3=29.65
1050 real, parameter :: svpt0=273.15
1051 real, parameter :: eps = 0.622
1053 real startlat, startlon, deltalat, deltalon
1055 call get_storage(iiplvl, 'P', P, ix, jx)
1056 call get_storage(iiplvl, 'TT', T, ix, jx)
1057 call get_storage(iiplvl, 'SH', Q, ix, jx)
1059 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1061 call put_storage(iiplvl, 'RH', rh, ix, jx)
1063 end subroutine g2_compute_rh_spechumd_upa
1065 !*****************************************************************************!
1067 SUBROUTINE g2_compute_rh_spechumd_upa2(ix, jx, iiplvl)
1068 ! Compute relative humidity from specific humidity in the upper air.
1069 use storage_module
1070 implicit none
1071 integer :: ix, jx
1072 integer :: iiplvl
1073 real :: lat1, lon1, dx, dy
1074 real, dimension(ix,jx) :: T, TH, RH, Q, P
1076 real, parameter :: svp1=611.2
1077 real, parameter :: svp2=17.67
1078 real, parameter :: svp3=29.65
1079 real, parameter :: svpt0=273.15
1080 real, parameter :: eps = 0.622
1082 real startlat, startlon, deltalat, deltalon
1084 call get_storage(iiplvl, 'THETA', TH, ix, jx)
1085 call get_storage(iiplvl, 'TT', T, ix, jx)
1086 call get_storage(iiplvl, 'SPECHUMD', Q, ix, jx)
1088 p=1.e5*(t/th)**(1005/287.05)
1090 rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
1092 call put_storage(iiplvl, 'RH', rh, ix, jx)
1094 end subroutine g2_compute_rh_spechumd_upa2
1096 !*****************************************************************************!
1098 SUBROUTINE g2_compute_pressure_tth_upa(ix, jx, iiplvl)
1099 ! Compute relative humidity from specific humidity in the upper air.
1100 use storage_module
1101 implicit none
1102 integer :: ix, jx
1103 integer :: iiplvl
1104 real :: lat1, lon1, dx, dy
1105 real, dimension(ix,jx) :: T, TH, P
1107 real, parameter :: svp1=611.2
1108 real, parameter :: svp2=17.67
1109 real, parameter :: svp3=29.65
1110 real, parameter :: svpt0=273.15
1111 real, parameter :: eps = 0.622
1113 real startlat, startlon, deltalat, deltalon
1115 call get_storage(iiplvl, 'THETA', TH, ix, jx)
1116 call get_storage(iiplvl, 'TT', T, ix, jx)
1118 p=1.e5*(t/th)**(1005/287.05)
1120 call put_storage(iiplvl, 'PRESSURE', p, ix, jx)
1122 end subroutine g2_compute_pressure_tth_upa
1124 !*****************************************************************************!
1126 subroutine ncep_grid_num (pnum)
1127 !
1128 ! Find the grib number for descriptive labelling.
1129 ! Grib2 doesn't have a grid-number entry, so we have to figure it out
1130 ! from the parameters.
1131 !
1132 use gridinfo ! Included to define map%
1133 integer :: pnum
1134 real, parameter :: eps = .01
1135 character (len=8) :: tmp8
1137 ! write(6,*) 'begin ncep_grid_num'
1138 ! write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
1139 tmp8 = ' '
1140 if (pnum .eq. 30) then ! lambert conformal
1141 if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
1142 write(tmp8,'("GRID 218")')
1143 else if (abs(map%dx - 40.63525) .lt. eps
1144 & .and. map%nx .eq. 185) then
1145 write(tmp8,'("GRID 212")')
1146 else if (abs(map%dx - 40.63525) .lt. eps
1147 & .and. map%nx .eq. 151) then
1148 write(tmp8,'("GRID 236")')
1149 else if (abs(map%dx - 81.2705) .lt. eps
1150 & .and. map%nx .eq. 93) then
1151 write(tmp8,'("GRID 211")')
1152 else if (abs (map%dx - 32.46341) .lt. eps
1153 & .and. map%nx .eq. 349) then
1154 write(tmp8,'("GRID 221")')
1155 else if (abs(map%dx - 20.317625) .lt. eps
1156 & .and. map%nx .eq. 301) then
1157 write(tmp8,'("GRID 252")')
1158 else if (abs(map%dx - 13.545087) .lt. eps
1159 & .and. map%nx .eq. 451) then
1160 write(tmp8,'("GRID 130")')
1161 endif
1162 else if (pnum .eq. 20) then ! polar stereographic
1163 if (abs(map%dx - 15.0) .lt. eps) then
1164 write(tmp8,'("GRID 88")')
1165 endif
1166 else if (pnum .eq. 0) then ! lat/lon
1167 if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
1168 write(tmp8,'("GRID 3")')
1169 else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
1170 write(tmp8,'("GRID 4")')
1171 endif
1172 endif
1173 map%source(25:32) = tmp8
1174 ! write(6,*) 'map%source = ',map%source
1175 end subroutine ncep_grid_num
1176 !*****************************************************************************!
1178 function earth_radius (icode, iscale, irad_m)
1179 ! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
1180 use module_debug
1181 real :: earth_radius
1182 integer :: icode
1183 integer :: iscale, irad_m
1184 if ( icode .eq. 0 ) then
1185 earth_radius = 6367470. * .001
1186 else if ( icode .eq. 1) then
1187 earth_radius = 0.001 * float(irad_m) / 10**iscale
1188 else if ( icode .eq. 6 ) then
1189 earth_radius = 6371229. * .001
1190 else if ( icode .eq. 8 ) then
1191 earth_radius = 6371200. * .001
1192 else
1193 call mprintf(.true.,ERROR,
1194 & "unknown earth radius for code %i",newline=.true.,i1=icode)
1195 endif
1196 end function earth_radius