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