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