ungrib/src/g2print.F

   1 !*****************************************************************************!
   2   program g2print                                                             !
   3 !                                                                             !
   4   use table
   5   use gridinfo
   6   use filelist
   7   implicit none
   8   interface
   9      subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2,&
  10           a3, h3, i3, l3, hlast)
  11        integer :: err
  12        character(len=*) , optional :: a1, a2, a3
  13        character(len=*), optional :: h1, h2, h3
  14        integer , optional :: i1, i2, i3
  15        logical, optional :: l1, l2, l3
  16        character(len=*), optional :: hlast
  17      end subroutine parse_args
  18   end interface
  19
  20   integer :: nunit1 = 12
  21   character(LEN=120) :: gribflnm
  22
  23   integer :: iprint
  24
  25   integer , parameter :: maxlvl = 100
  26
  27   real :: startlat, startlon, deltalat, deltalon
  28   real :: level
  29   character (LEN=9) ::  field
  30   character (LEN=3) ::  out_format
  31
  32   logical :: readit
  33
  34   integer, dimension(255) :: iuarr = 0
  35
  36   character (LEN=19) :: HSTART, HEND, HDATE
  37   character(LEN=19) :: hsave  = '0000-00-00_00:00:00'
  38   integer :: itime
  39   integer :: ntimes
  40   integer :: interval
  41   integer :: ierr
  42   logical :: ordered_by_date
  43   integer :: debug_level
  44   integer :: grib_version
  45   integer :: vtable_columns
  46
  47   character(len=30) :: hopt
  48   logical :: ivb = .FALSE.
  49   logical :: idb = .FALSE.
  50
  51 ! -----------------
  52
  53   gribflnm = '  '
  54   call parse_args(ierr, a1='v', l1=ivb, a2='V', l2=idb, hlast=gribflnm)
  55   if (ierr.ne.0) then
  56      call getarg(0, hopt)
  57      write(*,'(//,"Usage: ", A, " [-v] [-V] file",/)') trim(hopt)
  58      write(*,'("     -v   : Print more information about the GRIB records")')
  59      write(*,'("     -V   : Print way too much information about the GRIB&
  60           & records")')
  61      write(*,'("     file : GRIB file to read"//)')
  62      stop
  63   endif
  64
  65 ! -----------------
  66 ! Determine GRIB Edition number
  67   grib_version=0
  68   call edition_num(nunit1, trim(gribflnm), grib_version, ierr)
  69   if (ierr.eq.3) STOP 'GRIB file problem'
  70
  71
  72      debug_level = 0
  73      if (ivb) debug_level = 51
  74      if (idb) debug_level = 101
  75      write(6,*) 'reading from grib file = ',gribflnm
  76
  77      LOOP1 : DO
  78         ! At the beginning of LOOP1, we are at a new time period.
  79         ! Clear the storage arrays and associated level information.
  80
  81            ! If we need to read a new grib record, then read one.
  82
  83               if (grib_version.ne.2) then
  84 !                write(6,*) 'calling r_grib1 with iunit ', nunit1
  85 !                write(6,*) 'flnm = ',gribflnm
  86                  write(6,*) 'This is a Grib1 file. Please use g1print.\n'
  87                  stop
  88                  ! Read one record at a time from GRIB1 (and older Editions)
  89 !                call r_grib1(nunit1, gribflnm, level, field, &
  90 !                     hdate, debug_level, ierr, iuarr, iprint)
  91               else
  92
  93                  ! Read one file of records from GRIB2.
  94                  if (debug_level .gt. 100) write(6,*) 'calling r_grib2'
  95                  call r_grib2(nunit1, gribflnm, hdate, &
  96                       grib_version, debug_level, ierr)
  97
  98               endif
  99
 100               if (ierr.eq.1) then
 101                  ! We have hit the end of a file.  Exit LOOP1.
 102                  exit LOOP1
 103               endif
 104
 105      enddo LOOP1
 106
 107      if (grib_version.ne.2) then
 108         call c_close(iuarr(nunit1), iprint, ierr)
 109         iuarr(nunit1) = 0
 110      endif
 111
 112 ! And Now we are done:
 113
 114    print*,' '
 115    print*,' '
 116    print*,'  Successful completion of g2print   '
 117
 118 contains
 119   subroutine sort_filedates
 120     implicit none
 121
 122     integer :: n
 123     logical :: done
 124     if (nfiles > 1) then
 125        done = .FALSE.
 126        do while ( .not. done)
 127           done = .TRUE.
 128           do n = 1, nfiles-1
 129              if (filedates(n) > filedates(n+1)) then
 130                 filedates(size(filedates)) = filedates(n)
 131                 filedates(n) = filedates(n+1)
 132                 filedates(n+1) = filedates(size(filedates))
 133                 filedates(size(filedates)) = '0000-00-00 00:00:00.0000'
 134                 done = .FALSE.
 135              endif
 136           enddo
 137        enddo
 138     endif
 139   end subroutine sort_filedates
 140
 141 end program g2print
 142
 143 !*****************************************************************************!
 144
 145       SUBROUTINE r_grib2(junit, gribflnm, hdate,  &
 146         grib_edition, debug_level, ireaderr)
 147
 148       use grib_mod
 149       use params
 150       use table          ! Included to define g2code
 151       use gridinfo       ! Included to define map%
 152
 153       real, allocatable, dimension(:) :: hold_array
 154       parameter(msk1=32000,msk2=4000)
 155       character(len=1),allocatable,dimension(:) :: cgrib
 156       integer :: listsec0(3)
 157       integer :: listsec1(13)
 158       integer year, month, day, hour, minute, second, fcst
 159       character(len=*)  :: gribflnm
 160       character(len=*)  :: hdate
 161       character(len=8)  :: pabbrev
 162       character(len=20) :: labbrev
 163       character(len=80) :: tabbrev
 164       integer :: lskip, lgrib
 165       integer :: junit, itot, icount, iseek
 166       integer :: grib_edition
 167       integer :: i, j, ireaderr, ith
 168       integer :: currlen
 169       logical :: unpack, expand
 170       type(gribfield) :: gfld
 171       real :: level
 172       real :: scale_factor
 173       integer :: iplvl, lvl2
 174       ! For subroutine outout
 175       integer , parameter :: maxlvl = 100
 176       real , dimension(maxlvl) :: plvl
 177       integer :: nlvl
 178       integer , dimension(maxlvl) :: level_array
 179       logical :: verbose=.false.
 180       logical :: first = .true.
 181       integer :: debug_level
 182       character(len=4) :: tmp4
 183       character(len=40) :: string
 184       character(len=13) :: pstring = ',t50,":",i14)'
 185       character(len=15) :: rstring = ',t50,":",f14.5)'
 186       character(len=13) :: astring = ',t50,":",a14)'
 187
 188 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 189 !  SET ARGUMENTS
 190
 191       if (debug_level .gt. 50 ) then
 192         unpack=.true.
 193       else
 194         unpack=.false.
 195       endif
 196       expand=.true.
 197       hdate = '0000-00-00_00:00:00'
 198       ierr=0
 199       itot=0
 200       icount=0
 201       iseek=0
 202       lskip=0
 203       lgrib=0
 204       currlen=0
 205       ith=1
 206       scale_factor = 1e6
 207 !     do j = 1,10
 208 !       write(6,'("j = ",i4," level1 = ",i8," level2 = ",i8)') j, &
 209 !         level1(j), level2(j)
 210 !     enddo
 211
 212 !/* IOS Return Codes from BACIO:  */
 213 !/*  0    All was well                                   */
 214 !/* -1    Tried to open read only _and_ write only       */
 215 !/* -2    Tried to read and write in the same call       */
 216 !/* -3    Internal failure in name processing            */
 217 !/* -4    Failure in opening file                        */
 218 !/* -5    Tried to read on a write-only file             */
 219 !/* -6    Failed in read to find the 'start' location    */
 220 !/* -7    Tried to write to a read only file             */
 221 !/* -8    Failed in write to find the 'start' location   */
 222 !/* -9    Error in close                                 */
 223 !/* -10   Read or wrote fewer data than requested        */
 224
 225 !if ireaderr =1 we have hit the end of a file.
 226 !if ireaderr =2 we have hit the end of all the files.
 227
 228
 229       if ( debug_level .gt. 100 ) verbose = .true.
 230       if (verbose) write(6,*) 'begin r_grib2, flnm = ',gribflnm
 231       ! Open a byte-addressable file.
 232       CALL BAOPENR(junit,gribflnm,IOS)
 233       first = .true.
 234       if (verbose) write(6,*) 'back from baopenr, ios = ',ios
 235       if (ios.eq.0) then
 236       VERSION: do
 237
 238          ! Search opend file for the next GRIB2 messege (record).
 239       if (verbose) write(6,*) 'calling skgb'
 240          call skgb(junit,iseek,msk1,lskip,lgrib)
 241
 242          ! Check for EOF, or problem
 243          if (lgrib.eq.0) then
 244             exit
 245          endif
 246
 247          ! Check size, if needed allocate more memory.
 248          if (lgrib.gt.currlen) then
 249             if (allocated(cgrib)) deallocate(cgrib)
 250             allocate(cgrib(lgrib),stat=is)
 251             !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
 252             currlen=lgrib
 253          endif
 254
 255          ! Read a given number of bytes from unblocked file.
 256          call baread(junit,lskip,lgrib,lengrib,cgrib)
 257
 258          if (lgrib.ne.lengrib) then
 259             print *,'G2 r_grib2: IO Error.',lgrib,".ne.",lengrib
 260             stop 9
 261          endif
 262          iseek=lskip+lgrib
 263          icount=icount+1
 264
 265          if (verbose)  PRINT *,'G2 GRIB MESSAGE ',icount,' starts at',lskip+1
 266
 267          ! Unpack GRIB2 field
 268          call gb_info(cgrib,lengrib,listsec0,listsec1, &
 269                       numfields,numlocal,maxlocal,ierr)
 270          if (ierr.ne.0) then
 271            write(*,*) ' ERROR querying GRIB2 message = ',ierr
 272            stop 10
 273          endif
 274          itot=itot+numfields
 275
 276          grib_edition=listsec0(2)
 277          if (grib_edition.ne.2) then
 278               exit VERSION
 279          endif
 280
 281          ! Additional print statments for developer.
 282          if (verbose) then
 283            print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
 284            print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
 285            print *,'G2 Contains ',numlocal,' Local Sections ', &
 286                  ' and ',numfields,' data fields.'
 287          endif
 288
 289          ! ----
 290          ! Once per file fill in date, model and projection values.
 291
 292          if (first) then
 293            first = .false.
 294
 295            ! Build the 19-character date string, based on GRIB2 header date
 296            ! and time information, including forecast time information:
 297
 298            n=1
 299            call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
 300
 301           if (debug_level .gt. 100 ) then
 302           write(6,*) 'gfld%version = ',gfld%version
 303           if (gfld%discipline .eq. 0) then
 304             string = 'Meteorological products'
 305           else if (gfld%discipline .eq. 1) then
 306             string = 'Hydrological products'
 307           else if (gfld%discipline .eq. 2) then
 308             string = 'Land Surface products'
 309           else
 310             string = 'See code table 0.0'
 311           endif
 312           write(6,*) 'Discipline = ',gfld%discipline,'   ',string
 313           write(6,*) 'gfld%idsect(1) = ',gfld%idsect(1)
 314           write(6,*) 'gfld%idsect(2) = ',gfld%idsect(2)
 315           write(6,*) 'gfld%idsect(3) = ',gfld%idsect(3)
 316           write(6,*) 'gfld%idsect(4) = ',gfld%idsect(4)
 317           write(6,*) 'gfld%idsect(5) = ',gfld%idsect(5)
 318           write(6,*) 'gfld%idsect(6) = ',gfld%idsect(6)
 319           write(6,*) 'gfld%idsect(7) = ',gfld%idsect(7)
 320           write(6,*) 'gfld%idsect(8) = ',gfld%idsect(8)
 321           write(6,*) 'gfld%idsect(9) = ',gfld%idsect(9)
 322           write(6,*) 'gfld%idsect(10) = ',gfld%idsect(10)
 323           write(6,*) 'gfld%idsect(11) = ',gfld%idsect(11)
 324           write(6,*) 'gfld%idsect(12) = ',gfld%idsect(12)
 325           write(6,*) 'gfld%idsect(13) = ',gfld%idsect(13)
 326
 327           write(6,*) 'gfld%idsectlen = ',gfld%idsectlen
 328           write(6,*) 'gfld%locallen = ',gfld%locallen
 329           write(6,*) 'gfld%ifldnum = ',gfld%ifldnum
 330           write(6,*) 'gfld%ngrdpts = ',gfld%ngrdpts
 331           write(6,*) 'gfld%numoct_opt = ',gfld%numoct_opt
 332           write(6,*) 'gfld%interp_opt = ',gfld%interp_opt
 333
 334           write(6,*) 'gfld%griddef = ',gfld%griddef
 335           if (gfld%igdtnum .eq. 0) then
 336             string = 'Lat/Lon cylindrical equidistant'
 337           else if (gfld%igdtnum .eq. 1) then
 338             string = 'Rotated Lat/Lon'
 339           else if (gfld%igdtnum .eq. 2) then
 340             string = 'Stretched Lat/Lon'
 341           else if (gfld%igdtnum .eq. 20) then
 342             string = 'Polar Stereographic'
 343           else if (gfld%igdtnum .eq. 30) then
 344             string = 'Lambert Conformal'
 345           else if (gfld%igdtnum .eq. 40) then
 346             string = 'Gaussian Lat/Lon'
 347           else if (gfld%igdtnum .eq. 50) then
 348             string = 'Spherical harmonic coefficients'
 349           else
 350             string = 'see code table 3.1'
 351           endif
 352           write(6,*) 'Grid Template number = ',gfld%igdtnum,'   ',string
 353           write(6,*) 'gfld%igdtlen = ',gfld%igdtlen
 354           do i = 1, gfld%igdtlen
 355             write(6,*) 'gfld%igdtmpl(',i,') = ',gfld%igdtmpl(i)
 356           enddo
 357
 358           write(6,*) 'gfld%ipdtnum = ',gfld%ipdtnum
 359           write(6,*) 'gfld%ipdtlen = ',gfld%ipdtlen
 360           if ( gfld%ipdtnum .eq. 0 ) then
 361             do i = 1, gfld%ipdtlen
 362               write(6,*) 'gfld%ipdtmpl(',i,') = ',gfld%ipdtmpl(i)
 363             enddo
 364           endif
 365           write(6,*) 'gfld%num_coord = ',gfld%num_coord
 366           write(6,*) 'gfld%ndpts = ',gfld%ndpts
 367           write(6,*) 'gfld%idrtnum = ',gfld%idrtnum
 368           write(6,*) 'gfld%idrtlen = ',gfld%idrtlen
 369           write(6,*) 'gfld%expanded = ',gfld%expanded
 370           write(6,*) 'gfld%ibmap = ',gfld%ibmap
 371           endif
 372
 373           if (debug_level .le. 50) then
 374           write(6,*) '---------------------------------------------------------------------------------------'
 375           write(6,*) ' rec Prod Cat Param  Lvl    Lvl      Lvl     Prod    Name            Time          Fcst'
 376           write(6,*) ' num Disc     num    code   one      two     Templ                                 hour'
 377           write(6,*) '---------------------------------------------------------------------------------------'
 378           endif
 379
 380            year  =gfld%idsect(6)     !(FOUR-DIGIT) YEAR OF THE DATA
 381            month =gfld%idsect(7)     ! MONTH OF THE DATA
 382            day   =gfld%idsect(8)     ! DAY OF THE DATA
 383            hour  =gfld%idsect(9)     ! HOUR OF THE DATA
 384            minute=gfld%idsect(10)    ! MINUTE OF THE DATA
 385            second=gfld%idsect(11)    ! SECOND OF THE DATA
 386
 387            fcst = 0
 388
 389            ! Extract forecast time.
 390            if ( gfld%ipdtmpl(8) .eq. 1 ) then   ! time units are hours
 391              fcst = gfld%ipdtmpl(9)
 392            else if ( gfld%ipdtmpl(8) .eq. 0 ) then  ! minutes
 393              fcst = gfld%ipdtmpl(9) / 60.
 394            else if ( gfld%ipdtmpl(8) .eq. 2 ) then  ! days
 395              fcst = gfld%ipdtmpl(9) * 24.
 396            else
 397              fcst = 999
 398            endif
 399
 400            ! Compute valid time.
 401
 402           if (verbose) then
 403             print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
 404             print *, 'hhmm  ',gfld%idsect(9),gfld%idsect(10)
 405           endif
 406
 407            call build_hdate(hdate,year,month,day,hour,minute,second)
 408            if (verbose) print *, 'G2 hdate = ',hdate
 409 !          call geth_newdate(hdate,hdate,3600*fcst)   ! no need for thin in print
 410 !          print *, 'G2 hdate (fcst?) = ',hdate
 411
 412            !--
 413
 414            ! Indicator of the source (center) of the data.
 415            icenter = gfld%idsect(1)
 416
 417            ! Indicator of model (or whatever) which generated the data.
 418            iprocess = gfld%ipdtmpl(5)
 419
 420
 421            if (icenter.eq.7) then
 422              if (iprocess.eq.83 .or. iprocess.eq.84) then
 423                map%source = 'NCEP NAM Model'
 424              elseif (iprocess.eq.81) then
 425                map%source = 'NCEP GFS Model'
 426              elseif (iprocess.eq.82) then
 427                map%source = 'NCEP GFS GDAS'
 428              elseif (iprocess.eq.96) then
 429                map%source = 'NCEP GFS Model'
 430              elseif (iprocess.eq.86 .or. iprocess.eq.100) then
 431                map%source = 'NCEP RUC Model'    ! 60 km
 432              elseif (iprocess.eq.101) then
 433                map%source = 'NCEP RUC Model'    ! 40 km
 434              elseif (iprocess.eq.105) then
 435                map%source = 'NCEP RUC Model'    ! 20 km
 436              elseif (iprocess.eq.107) then
 437                map%source = 'NCEP GEFS'
 438              elseif (iprocess.eq.109) then
 439                map%source = 'NCEP RTMA'
 440              elseif (iprocess.eq.140) then
 441                map%source = 'NCEP NARR'
 442              elseif (iprocess.eq.44) then
 443                map%source = 'NCEP SST Analysis'
 444              elseif (iprocess.eq.70) then
 445                map%source = 'GFDL Hurricane Model'
 446              elseif (iprocess.eq.80) then
 447                map%source = 'NCEP GFS Ensemble'
 448              elseif (iprocess.eq.107) then          ! renumbered as of 23 Feb 2010
 449                map%source = 'NCEP GFS Ensemble'
 450              elseif (iprocess.eq.129) then
 451                map%source = 'NCEP GODAS'
 452              elseif (iprocess.eq.25) then
 453                map%source = 'NCEP SNOW COVER ANALYSIS'
 454              else
 455                map%source = 'unknown model from NCEP'
 456                write (6,*) 'unknown NCEP model, iprocess = ',iprocess
 457              end if
 458            else if (icenter .eq. 57) then
 459              if (iprocess .eq. 87) then
 460                map%source = 'AFWA AGRMET'
 461              else
 462                map%source = 'AFWA'
 463              endif
 464            else if (icenter .eq. 58) then
 465                map%source = 'US Navy FNOC'
 466            else if (icenter .eq. 98) then
 467              map%source = 'ECMWF'
 468            else if (icenter .eq. 34) then
 469              map%source = 'JMA'
 470            else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
 471              map%source = 'UKMO'
 472            else
 473              map%source = 'unknown model and orig center'
 474            end if
 475
 476            !--
 477
 478            ! Store information about the grid on which the data is.
 479            ! This stuff gets stored in the MAP variable, as defined in
 480            ! module GRIDINFO.
 481
 482            map%startloc = 'SWCORNER'
 483
 484            if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
 485               map%igrid = 0
 486               map%nx = gfld%igdtmpl(8)
 487               map%ny = gfld%igdtmpl(9)
 488               map%dx = gfld%igdtmpl(17)
 489               map%dy = gfld%igdtmpl(18)
 490               map%lat1 = gfld%igdtmpl(12)
 491               map%lon1 = gfld%igdtmpl(13)
 492
 493               if ((gfld%igdtmpl(10) .eq. 0).OR.   &
 494                   (gfld%igdtmpl(10) .eq. 255)) THEN
 495           ! Scale lat/lon values to 0-180, default range is 1e6.
 496                 map%lat1 = map%lat1/scale_factor
 497                 map%lon1 = map%lon1/scale_factor
 498           ! Scale dx/dy values to degrees, default range is 1e6.
 499                 map%dx = map%dx/scale_factor
 500                 map%dy = map%dy/scale_factor
 501               else
 502           ! Basic angle and subdivisions are non-zero (not tested)
 503                 map%lat1 = map%lat1 *  &
 504                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 505                 map%lon1 = map%lon1 *  &
 506                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 507                 map%dx = map%dx *  &
 508                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 509                 map%dy = map%dy *  &
 510                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 511               endif
 512
 513            elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
 514               map%igrid = 3
 515               map%nx = gfld%igdtmpl(8)
 516               map%ny = gfld%igdtmpl(9)
 517               map%lov = gfld%igdtmpl(14)
 518               map%truelat1 = gfld%igdtmpl(19)
 519               map%truelat2 = gfld%igdtmpl(20)
 520               map%dx = gfld%igdtmpl(15)
 521               map%dy = gfld%igdtmpl(16)
 522               map%lat1 = gfld%igdtmpl(10)
 523               map%lon1 = gfld%igdtmpl(11)
 524
 525            elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
 526               map%igrid = 0
 527               map%nx = gfld%igdtmpl(8)
 528               map%ny = gfld%igdtmpl(9)
 529               map%dx = gfld%igdtmpl(17)
 530               map%dy = gfld%igdtmpl(18) ! ?not in Grid Definition Template 3.40 doc
 531               map%lat1 = gfld%igdtmpl(12)
 532               map%lon1 = gfld%igdtmpl(13)
 533
 534               ! Scale dx/dy values to degrees, default range is 1e6.
 535               if (map%dx.gt.10000) then
 536                  map%dx = map%dx/scale_factor
 537               endif
 538               if (map%dy.gt.10000) then
 539                  map%dy = (map%dy/scale_factor)*(-1)
 540               endif
 541
 542               ! Scale lat/lon values to 0-180, default range is 1e6.
 543               if (map%lat1.ge.scale_factor) then
 544                  map%lat1 = map%lat1/scale_factor
 545               endif
 546               if (map%lon1.ge.scale_factor) then
 547                  map%lon1 = map%lon1/scale_factor
 548               endif
 549            print *,'Gaussian Grid: Dx,Dy,lat,lon',map%dx,map%dy, &
 550              map%lat1,map%lon1
 551
 552            elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
 553               map%igrid = 5
 554               map%nx = gfld%igdtmpl(8)
 555               map%ny = gfld%igdtmpl(9)
 556               !map%lov = gfld%igdtmpl(x) ! ?not in Grid Definition Template 3.20 doc
 557               map%truelat1 = 60.
 558               map%truelat2 = 91.
 559               !map%dx = gfld%igdtmpl(x)
 560               !map%dy = gfld%igdtmpl(x)
 561               map%lat1 = gfld%igdtmpl(10)
 562               map%lon1 = gfld%igdtmpl(11)
 563
 564            else
 565               print*, 'GRIB2 Unknown Projection: ',gfld%igdtnum
 566               print*, 'see Code Table 3.1: Grid Definition Template No'
 567            endif
 568
 569             call gf_free(gfld)
 570          endif
 571
 572          ! ----
 573
 574          ! Continue to unpack GRIB2 field.
 575          if (debug_level .gt. 100) write(6,*) 'numfields = ',numfields
 576          do n=1,numfields ! e.g. U and V would =2, otherwise its usually =1
 577            call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
 578            if (ierr.ne.0) then
 579              write(*,*) ' ERROR extracting field gf_getfld = ',ierr
 580              cycle
 581            endif
 582
 583 ! ------------------------------------
 584          ! Additional print information for developer.
 585          pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),  &
 586                                   gfld%ipdtmpl(2))
 587          if (debug_level .gt. 50 ) then
 588            print *
 589 !          print *,'G2 FIELD ',n
 590            if (n==1) then
 591             write(*,'(/,"GRIB2 SECTION 0 - INDICATOR SECTION:")')
 592             write(*,'(5x,"Discipline"'//pstring) gfld%discipline
 593             write(*,'(5x,"GRIB Edition Number"'//pstring) gfld%version
 594             write(*,'(5x,"GRIB length"'//pstring) lengrib
 595             write(*,'(/,"GRIB2 SECTION 1 - IDENTIFICATION SECTION:")')
 596             write(*,'(5x,"Length of Section"'//pstring) gfld%idsectlen
 597             write(*,'(5x,"Originating Center ID"'//pstring) &
 598                gfld%idsect(1)
 599             write(*,'(5x,"Subcenter ID"'//pstring) gfld%idsect(2)
 600             write(*,'(5x,"GRIB Master Table Version"'//pstring) &
 601             gfld%idsect(3)
 602             write(*,'(5x,"GRIB Local Table Version"'//pstring) &
 603             gfld%idsect(4)
 604             write(*,'(5x,"Significance of Reference Time"'//pstring) &
 605             gfld%idsect(5)
 606             write(*,'(5x,"Year"'//pstring)  gfld%idsect(6)
 607             write(*,'(5x,"Month"'//pstring)  gfld%idsect(7)
 608             write(*,'(5x,"Day"'//pstring)  gfld%idsect(8)
 609             write(*,'(5x,"Hour"'//pstring)  gfld%idsect(9)
 610             write(*,'(5x,"Minute"'//pstring)  gfld%idsect(10)
 611             write(*,'(5x,"Second"'//pstring)  gfld%idsect(11)
 612             write(*,'(5x,"Production Status of data"'//pstring) &
 613               gfld%idsect(12)
 614             write(*,'(5x,"Type of processed data"'//pstring) &
 615               gfld%idsect(13)
 616 !           print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
 617            endif
 618            write(*,'(/,"GRIB2 SECTION 2 - LOCAL SECTION:")')
 619            write(*,'(5x,"Length of Section 2"'//pstring) gfld%locallen
 620            if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
 621            do j = 1, gfld%locallen
 622              write(*,'(5x,"Local value "'//astring) gfld%local(j)
 623            enddo
 624 !             print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
 625            endif
 626            write(*,'(/,"GRIB2 SECTION 3 - GRID DEFINITION SECTION:")')
 627 !          write(*,'(5x,"Length of Section 3"'//pstring) gfld%unknown
 628            write(*,'(5x,"Source of grid definition"'&
 629            //pstring) gfld%griddef
 630            write(*,'(5x,"Number of grid points"'//pstring) gfld%ngrdpts
 631            write(*,'(5x,"Number of octets for addnl points"'//pstring) &
 632            gfld%numoct_opt
 633            write(*,'(5x,"Interpretation list"'//pstring) &
 634            gfld%interp_opt
 635            write(*,'(5x,"Grid Definition Template Number"'//pstring) &
 636            gfld%igdtnum
 637            if (gfld%igdtnum .eq. 0 .or. gfld%igdtnum .eq. 1 .or.  &
 638                gfld%igdtnum .eq. 2 .or. gfld%igdtnum .eq. 3 ) then
 639              if (gfld%igdtnum .eq. 0 ) then
 640                write(*,'(5x,"Lat/Lon or Cylindrical Equidistant Grid")')
 641              else if (gfld%igdtnum .eq. 1 ) then
 642                write(*,'(5x,"Rotated Lat/Lon or Cylind. Equi. Grid")')
 643              else if (gfld%igdtnum .eq. 2 ) then
 644                write(*,'(5x,"Stretched Lat/Lon or Cylind. Equi. Grid")')
 645              else if (gfld%igdtnum .eq. 3 ) then
 646                write(*,'(5x,"Stretched and Rotated Lat/Lon Grid")')
 647              endif
 648              write(*,'(10x,"Shape of the Earth"'//pstring) &
 649                 gfld%igdtmpl(1)
 650              write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
 651                 gfld%igdtmpl(2)
 652              write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
 653                 gfld%igdtmpl(3)
 654              write(*,'(10x,"Scale factor of major axis"'//pstring) &
 655                 gfld%igdtmpl(4)
 656              write(*,'(10x,"Scaled value of major axis"'//pstring) &
 657                 gfld%igdtmpl(5)
 658              write(*,'(10x,"Scale factor of minor axis"'//pstring) &
 659                 gfld%igdtmpl(6)
 660              write(*,'(10x,"Scaled value of minor axis"'//pstring) &
 661                 gfld%igdtmpl(7)
 662              write(*,'(10x,"Ni - points along a parallel"'//pstring) &
 663                  gfld%igdtmpl(8)
 664              write(*,'(10x,"Nj - points along a meridian"'//pstring) &
 665                  gfld%igdtmpl(9)
 666              write(*,'(10x,"Basic angle of initial domain"'//pstring)&
 667                  gfld%igdtmpl(10)
 668              write(*,'(10x,"Subdivisions of basic angle"'//pstring) &
 669                  gfld%igdtmpl(11)
 670              write(*,'(10x,"La1"'//pstring)  gfld%igdtmpl(12)
 671              write(*,'(10x,"Lo1"'//pstring)  gfld%igdtmpl(13)
 672              write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
 673               gfld%igdtmpl(14)
 674              write(*,'(10x,"La2"'//pstring) gfld%igdtmpl(15)
 675              write(*,'(10x,"Lo2"'//pstring) gfld%igdtmpl(16)
 676              write(*,'(10x,"Di - i-dir increment"'//pstring) &
 677                gfld%igdtmpl(17)
 678              write(*,'(10x,"Dj - j-dir increment"'//pstring) &
 679                gfld%igdtmpl(18)
 680              write(*,'(10x,"Scanning mode"'//pstring) &
 681                gfld%igdtmpl(19)
 682              if (gfld%igdtnum .eq. 1) then
 683              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 684                gfld%igdtmpl(20)
 685              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 686                gfld%igdtmpl(21)
 687              write(*,'(10x,"Angle of rotation of projection"'//pstring)&
 688                gfld%igdtmpl(22)
 689              else if (gfld%igdtnum .eq. 2) then
 690              write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
 691                gfld%igdtmpl(20)
 692              write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
 693                gfld%igdtmpl(21)
 694              write(*,'(10x,"Stretching factor"'//pstring) &
 695                gfld%igdtmpl(22)
 696              else if (gfld%igdtnum .eq. 3) then
 697              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 698                gfld%igdtmpl(20)
 699              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 700                gfld%igdtmpl(21)
 701              write(*,'(10x,"Angle of rotation of projection"'//pstring)&
 702                gfld%igdtmpl(22)
 703              write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
 704                gfld%igdtmpl(23)
 705              write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
 706                gfld%igdtmpl(24)
 707              write(*,'(10x,"Stretching factor"'//pstring) &
 708                gfld%igdtmpl(25)
 709              endif
 710            else if (gfld%igdtnum .eq. 10) then
 711              write(*,'(5x,"Mercator Grid")')
 712            else if (gfld%igdtnum .eq. 20 .or. gfld%igdtnum .eq. 30) then
 713              if (gfld%igdtnum .eq. 20) then
 714                write(*,'(5x,"Polar Stereographic Grid")')
 715              else if (gfld%igdtnum .eq. 30) then
 716                write(*,'(5x,"Lambert Conformal Grid")')
 717              endif
 718              write(*,'(10x,"Shape of the Earth"'//pstring) &
 719                 gfld%igdtmpl(1)
 720              write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
 721                 gfld%igdtmpl(2)
 722              write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
 723                 gfld%igdtmpl(3)
 724              write(*,'(10x,"Scale factor of major axis"'//pstring) &
 725                 gfld%igdtmpl(4)
 726              write(*,'(10x,"Scaled value of major axis"'//pstring) &
 727                 gfld%igdtmpl(5)
 728              write(*,'(10x,"Scale factor of minor axis"'//pstring) &
 729                 gfld%igdtmpl(6)
 730              write(*,'(10x,"Scaled value of minor axis"'//pstring) &
 731                 gfld%igdtmpl(7)
 732              write(*,'(10x,"Nx"'//pstring) gfld%igdtmpl(8)
 733              write(*,'(10x,"Ny"'//pstring) gfld%igdtmpl(9)
 734              write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(10)
 735              write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(11)
 736              write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
 737               gfld%igdtmpl(12)
 738              write(*,'(10x,"LaD"'//pstring) gfld%igdtmpl(13)
 739              write(*,'(10x,"LoV"'//pstring) gfld%igdtmpl(14)
 740              write(*,'(10x,"Dx"'//pstring) gfld%igdtmpl(15)
 741              write(*,'(10x,"Dy"'//pstring) gfld%igdtmpl(16)
 742              write(*,'(10x,"Projection Center Flag"'//pstring) &
 743                gfld%igdtmpl(17)
 744              write(*,'(10x,"Scanning mode"'//pstring) &
 745                gfld%igdtmpl(18)
 746              if (gfld%igdtnum .eq. 30) then
 747              write(*,'(10x,"Latin 1 "'//pstring) &
 748                gfld%igdtmpl(19)
 749              write(*,'(10x,"Latin 2 "'//pstring) &
 750                gfld%igdtmpl(20)
 751              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 752                gfld%igdtmpl(21)
 753              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 754                gfld%igdtmpl(22)
 755              endif
 756            else if (gfld%igdtnum .eq. 40 .or. gfld%igdtnum .eq. 41) then
 757              if (gfld%igdtnum .eq. 40) then
 758                write(*,'(5x,"Gaussian Lat/Lon Grid")')
 759              else if (gfld%igdtnum .eq. 41) then
 760                write(*,'(5x,"Rotated Gaussian Lat/Lon Grid")')
 761              else if (gfld%igdtnum .eq. 42) then
 762                write(*,'(5x,"Stretched Gaussian Lat/Lon Grid")')
 763              else if (gfld%igdtnum .eq. 43) then
 764                write(*,'(5x,"Stretched and Rotated Gaussian Lat/Lon ")')
 765              endif
 766            else
 767            do j = 1, gfld%igdtlen
 768              write(*,'(5x,"Grid Definition Template entry "'//pstring) &
 769              gfld%igdtmpl(j)
 770            enddo
 771            endif
 772 !          print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts, &
 773 !                                 gfld%numoct_opt,gfld%interp_opt,  &
 774 !                                 gfld%igdtnum
 775 !          print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',   &
 776 !                 (gfld%igdtmpl(j),j=1,gfld%igdtlen)
 777            if ( gfld%num_opt .eq. 0 ) then
 778 !            print *,'G2 NO Section 3 List Defining No. of Data Points.'
 779            else
 780              print *,'G2 Section 3 Optional List: ',     &
 781                       (gfld%list_opt(j),j=1,gfld%num_opt)
 782            endif
 783           write(*,'(/,"GRIB2 SECTION 4 - PRODUCT DEFINITION SECTION:")')
 784 !          write(*,'(5x,"Length of Section 4"'//pstring) gfld%unknown
 785            write(*,'(5x,"Product Definition Template Number"'//pstring)&
 786             gfld%ipdtnum
 787             do j = 1, gfld%ipdtlen
 788               write(tmp4,'(i4)') j
 789               string = '(5x,"Template Entry  '//tmp4 // '"'
 790               write(*,string//pstring) gfld%ipdtmpl(j)
 791             enddo
 792 !          print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',   &
 793 !               (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
 794
 795            !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
 796            !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
 797            write(*,'(5x,"Product Abbreviated Name",t50,":",a14)')&
 798            pabbrev
 799 !           print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
 800
 801            if ( gfld%num_coord .eq. 0 ) then
 802 !            print *,'G2 NO Optional Vertical Coordinate List.'
 803            else
 804              print *,'G2 Section 4 Optional Coordinates: ',   &
 805                    (gfld%coord_list(j),j=1,gfld%num_coord)
 806            endif
 807 !          if ( gfld%ibmap .ne. 255 ) then
 808 !             print *,'G2 Num. of Data Points = ',gfld%ndpts,   &
 809 !                  '    with BIT-MAP ',gfld%ibmap
 810 !          else
 811 !             print *,'G2 Num. of Data Points = ',gfld%ndpts,  &
 812 !                     '    NO BIT-MAP '
 813 !          endif
 814          write(*,'(/,"GRIB2 SECTION 5 - DATA REPRESENTATION SECTION:")')
 815          write(*,'(5x,"Data Representation Template Number"'//pstring)&
 816             gfld%idrtnum
 817             do j = 1, gfld%idrtlen
 818               write(tmp4,'(i4)') j
 819               string = '(5x,"Template Entry  '//tmp4 // '"'
 820               write(*,string//pstring) gfld%idrtmpl(j)
 821             enddo
 822 !          print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',  &
 823 !               (gfld%idrtmpl(j),j=1,gfld%idrtlen)
 824
 825 !      if ( gfld%ipdtnum .eq. 0 ) then
 826 !        if (gfld%ipdtmpl(1) .eq. 0 ) then
 827 !          write(6,*) 'Temperature'
 828 !        else if (gfld%ipdtmpl(1) .eq. 1 ) then
 829 !          write(6,*) 'Moisture'
 830 !        else if (gfld%ipdtmpl(1) .eq. 2 ) then
 831 !          write(6,*) 'Momentum'
 832 !        else if (gfld%ipdtmpl(1) .eq. 3 ) then
 833 !         write(6,*) 'Mass'
 834 !       endif
 835 !      endif
 836
 837          write(*,'(/,"GRIB2 SECTION 6 - BIT-MAP SECTION:")')
 838          write(*,'(5x,"Bit-map indicator"'//pstring) &
 839            gfld%ibmap
 840
 841            fldmax=gfld%fld(1)
 842            fldmin=gfld%fld(1)
 843            sum=gfld%fld(1)
 844            do j=2,gfld%ndpts
 845              if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
 846              if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
 847              sum=sum+gfld%fld(j)
 848            enddo ! gfld%ndpts
 849
 850          write(*,'(/,"GRIB2 SECTION 7 - DATA SECTION:")')
 851
 852          write(*,'(5x,"Minimum Data Value "'//rstring)&
 853             fldmin
 854          write(*,'(5x,"Maximum Data Value "'//rstring)&
 855             fldmax
 856 !          print *,'G2 Data Values:'
 857 !          write(*,fmt='("G2 MIN=",f21.8," AVE=",f21.8,    &
 858 !               " MAX=",f21.8)') fldmin,sum/gfld%ndpts,fldmax
 859          if (debug_level .gt. 100 ) then
 860             print*, 'gfld%fld = ',gfld%fld
 861 !           do j=1,gfld%ndpts
 862 !              write(*,*) j, gfld%fld(j)
 863 !           enddo
 864          endif
 865          endif ! Additional Print information
 866 ! ------------------------------------
 867          if ( debug_level .le. 50) then
 868            if(gfld%ipdtmpl(10).eq.100) then    ! pressure level
 869              level=gfld%ipdtmpl(12) *  &
 870                      (10. ** (-1. * gfld%ipdtmpl(11)))
 871            else if(gfld%ipdtmpl(10).eq.101 .or.&    ! sea level, sfc, or trop
 872                 gfld%ipdtmpl(10).eq.1 .or. gfld%ipdtmpl(10).eq.7) then
 873              level = 0
 874            else
 875              level=gfld%ipdtmpl(12)
 876            endif
 877            if (gfld%ipdtmpl(13) .eq. 255) then
 878              lvl2 = 0
 879            else
 880              lvl2 = gfld%ipdtmpl(15)
 881            endif
 882 !           Account for multiple forecast hours in one file
 883            if (gfld%ipdtnum.eq.0 .or. gfld%ipdtnum.eq.1 .or. gfld%ipdtnum.eq. 8) then
 884               ! Product Definition Template 4.0, 4.1, 4.8
 885               ! Extract forecast time.
 886              if ( gfld%ipdtmpl(8) .eq. 1 ) then   ! time units are hours
 887                fcst = gfld%ipdtmpl(9)
 888              else if ( gfld%ipdtmpl(8) .eq. 0 ) then  ! minutes
 889                fcst = gfld%ipdtmpl(9) / 60.
 890              else if ( gfld%ipdtmpl(8) .eq. 2 ) then  ! days
 891                fcst = gfld%ipdtmpl(9) * 24.
 892              else
 893                fcst = 999
 894              endif
 895            endif
 896
 897          write(6,987) itot,gfld%discipline,gfld%ipdtmpl(1),  &
 898            gfld%ipdtmpl(2),gfld%ipdtmpl(10),int(level),&
 899            lvl2,gfld%ipdtnum,pabbrev,hdate,fcst
 900   987     format(2i4,i5,i4,i8,i8,i8,i8,3x,a10,a20,i5.2)
 901          endif
 902
 903          ! Deallocate arrays decoding GRIB2 record.
 904          call gf_free(gfld)
 905
 906          enddo ! 1,numfields
 907
 908       enddo VERSION ! skgb
 909
 910        if (debug_level .gt. 50) &
 911           print *, 'G2 total number of fields found = ',itot
 912
 913         CALL BACLOSE(junit,IOS)
 914
 915         ireaderr=1
 916       else
 917         print *,'open status failed because',ios
 918         hdate = '9999-99-99_99:99:99'
 919         ireaderr=2
 920       endif ! ireaderr check
 921
 922       END subroutine r_grib2
 923
 924 !*****************************************************************************!
 925 ! Subroutine edition_num                                                         !
 926 !                                                                             !
 927 ! Purpose:                                                                    !
 928 !    Read one record from the input GRIB2 file.  Based on the information in  !
 929 !    the GRIB2 header and the user-defined Vtable, decide whether the field in!
 930 !    the GRIB2 record is one to process or to skip.  If the field is one we   !
 931 !    want to keep, extract the data from the GRIB2 record, and pass the data  !
 932 !    back to the calling routine.                                             !
 933 !                                                                             !
 934 ! Argument list:                                                              !
 935 !    Input:                                                                   !
 936 !       JUNIT   : "Unit Number" to open and read from.  Not really a Fortran  !
 937 !                 unit number, since we do not do Fortran I/O for the GRIB2   !
 938 !                 files.  Nor is it a UNIX File Descriptor returned from a C  !
 939 !                 OPEN statement.  It is really just an array index to the    !
 940 !                 array (IUARR) where the UNIX File Descriptor values are     !
 941 !                 stored.                                                     !
 942 !       GRIB2FILE: File name to open, if it is not already open.              !
 943 !                                                                             !
 944 !    Output:                                                                  !
 945 !       GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2                !
 946 !       IERR     : Error flag: 0 - no error on read from GRIB2 file.          !
 947 !                              1 - Hit the end of the GRIB2 file.             !
 948 !                              2 - The file GRIBFLNM we tried to open does    !
 949 !                                  not exist.                                 !
 950 ! Author: Paula McCaslin                                                      !
 951 ! NOAA/FSL                                                                    !
 952 ! Sept 2004                                                                   !
 953 !*****************************************************************************!
 954
 955       SUBROUTINE edition_num(junit, gribflnm, grib_edition, ireaderr)
 956
 957       use grib_mod
 958       use params
 959
 960       parameter(msk1=32000,msk2=4000)
 961       character(len=1),allocatable,dimension(:) :: cgrib
 962       integer :: listsec0(3)
 963       integer :: listsec1(13)
 964       character(len=*)  :: gribflnm
 965       integer :: lskip, lgrib
 966       integer :: junit
 967       integer :: grib_edition
 968       integer :: i, j, ireaderr
 969       integer :: currlen
 970
 971       character(len=4) :: ctemp
 972       character(len=4),parameter :: grib='GRIB',c7777='7777'
 973
 974 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 975 !  SET ARGUMENTS
 976
 977       itot=0
 978       icount=0
 979       iseek=0
 980       lskip=0
 981       lgrib=0
 982       currlen=0
 983
 984 !/* IOS Return Codes from BACIO:  */
 985 !/*  0    All was well                                   */
 986 !/* -1    Tried to open read only _and_ write only       */
 987 !/* -2    Tried to read and write in the same call       */
 988 !/* -3    Internal failure in name processing            */
 989 !/* -4    Failure in opening file                        */
 990 !/* -5    Tried to read on a write-only file             */
 991 !/* -6    Failed in read to find the 'start' location    */
 992 !/* -7    Tried to write to a read only file             */
 993 !/* -8    Failed in write to find the 'start' location   */
 994 !/* -9    Error in close                                 */
 995 !/* -10   Read or wrote fewer data than requested        */
 996
 997 !if ireaderr =1 we have hit the end of a file.
 998 !if ireaderr =2 we have hit the end of all the files.
 999 !if ireaderr =3 beginning characters 'GRIB' not found
1000
1001 !     write(6,*) 'junit = ',junit,' gribflnm = ',gribflnm
1002
1003       ! Open a byte-addressable file.
1004       CALL BAOPENR(junit,gribflnm,IOS)      ! from w3lib
1005 !     write(6,*) 'ios = ',ios
1006       if (ios.eq.0) then
1007
1008          ! Search opend file for the next GRIB2 messege (record).
1009          call skgb(junit,iseek,msk1,lskip,lgrib)
1010
1011          ! Check for EOF, or problem
1012          if (lgrib.eq.0) then
1013             write(*,'("\n\tThere is a problem with the input file.")')
1014             write(*,'("\tPerhaps it is not a Grib2 file?\n")')
1015             STOP "Grib2 file or date problem, stopping in edition_num."
1016          endif
1017
1018          ! Check size, if needed allocate more memory.
1019          if (lgrib.gt.currlen) then
1020             if (allocated(cgrib)) deallocate(cgrib)
1021             allocate(cgrib(lgrib),stat=is)
1022             currlen=lgrib
1023          endif
1024
1025          ! Read a given number of bytes from unblocked file.
1026          call baread(junit,lskip,lgrib,lengrib,cgrib)
1027
1028          ! Check for beginning of GRIB message in the first 100 bytes
1029          istart=0
1030          do j=1,100
1031             ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
1032             if (ctemp.eq.grib ) then
1033               istart=j
1034               exit
1035             endif
1036          enddo
1037          if (istart.eq.0) then
1038             ireaderr=3
1039             print*, "The beginning 4 characters >GRIB< were not found."
1040          endif
1041
1042          ! Unpack Section 0 - Indicator Section to extract GRIB edition field
1043          iofst=8*(istart+5)
1044          call gbyte(cgrib,discipline,iofst,8)     ! Discipline
1045          iofst=iofst+8
1046          call gbyte(cgrib,grib_edition,iofst,8)   ! GRIB edition number
1047
1048          print *, 'ungrib - grib edition num',  grib_edition
1049          CALL BACLOSE(junit,IOS)
1050          ireaderr=1
1051       else if (ios .eq. -4) then
1052         print *,'edition_num: unable to open ',gribflnm
1053         stop 'edition_num'
1054       else
1055          print *,'edition_num: open status failed because',ios,gribflnm
1056          ireaderr=2
1057       endif ! ireaderr check
1058
1059       END subroutine edition_num
1060 subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2, a3, h3, i3, l3, &
1061      hlast)
1062   integer :: err
1063   character(len=*) , optional :: a1, a2, a3
1064   character(len=*), optional :: h1, h2, h3
1065   integer , optional :: i1, i2, i3
1066   logical, optional :: l1, l2, l3
1067   character(len=*), optional :: hlast
1068
1069   character(len=100) :: hold
1070   integer :: ioff = 0
1071
1072   if (present(hlast)) then
1073      ioff = -1
1074   endif
1075
1076   err = 0
1077
1078   narg = iargc()
1079   numarg = narg + ioff
1080
1081   i = 1
1082   LOOP : do while ( i <= numarg)
1083
1084      ierr = 1
1085      if (present(i1)) then
1086         call checkiarg(i, a1, i1, ierr)
1087      elseif (present(h1)) then
1088         call checkharg(i, a1, h1, ierr)
1089      elseif (present(l1)) then
1090         call checklarg(i, a1, l1, ierr)
1091      endif
1092      if (ierr.eq.0) cycle LOOP
1093
1094      if (present(i2)) then
1095         call checkiarg(i, a2, i2, ierr)
1096      elseif (present(h2)) then
1097         call checkharg(i, a2, h2, ierr)
1098      elseif (present(l2)) then
1099         call checklarg(i, a2, l2, ierr)
1100      endif
1101      if (ierr.eq.0) cycle LOOP
1102
1103      if (present(i3)) then
1104         call checkiarg(i, a3, i3, ierr)
1105      elseif (present(h3)) then
1106         call checkharg(i, a3, h3, ierr)
1107      elseif (present(l3)) then
1108         call checklarg(i, a3, l3, ierr)
1109      endif
1110      if (ierr.eq.0) cycle LOOP
1111
1112      err = 1
1113      call getarg(1, hold)
1114      write(*, '("arg = ", A)') trim(hold)
1115
1116      exit LOOP
1117
1118   enddo LOOP
1119
1120   if (present(hlast)) then
1121      if (narg.eq.0) then
1122         err = 1
1123      else
1124         call getarg(narg, hlast)
1125      endif
1126   endif
1127
1128 contains
1129   subroutine checkiarg(c, a, i, ierr)
1130     integer :: c
1131     character(len=*) :: a
1132     integer :: i
1133
1134     character(len=100) :: hold
1135     ierr = 1
1136
1137     call getarg(c, hold)
1138
1139     if ('-'//a.eq.trim(hold)) then
1140        c = c + 1
1141        call getarg(c, hold)
1142        read(hold, *) i
1143        c = c + 1
1144        ierr = 0
1145     elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1146        hold = hold(len_trim(a)+2: len(hold))
1147        read(hold, *) i
1148        c = c + 1
1149        ierr = 0
1150     endif
1151
1152   end subroutine checkiarg
1153   subroutine checkharg(c, a, h, ierr)
1154     integer :: c
1155     character(len=*) :: a
1156     character(len=*) :: h
1157
1158     character(len=100) :: hold
1159     ierr = 1
1160
1161     call getarg(c, hold)
1162
1163     if ('-'//a.eq.trim(hold)) then
1164        c = c + 1
1165        call getarg(c, hold)
1166        h = trim(hold)
1167        c = c + 1
1168        ierr = 0
1169     elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1170        hold = hold(len_trim(a)+2: len(hold))
1171        h = trim(hold)
1172        c = c + 1
1173        ierr = 0
1174     endif
1175
1176   end subroutine checkharg
1177
1178   subroutine checklarg(c, a, l, ierr)
1179     integer :: c
1180     character(len=*) :: a
1181     logical :: l
1182
1183     character(len=100) :: hold
1184     ierr = 1
1185
1186     call getarg(c, hold)
1187     if ('-'//a.eq.trim(hold)) then
1188        l = .TRUE.
1189        c = c + 1
1190        ierr = 0
1191     endif
1192
1193   end subroutine checklarg
1194
1195 end subroutine parse_args
1196