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 = 150
  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 output
 175       integer , parameter :: maxlvl = 150
 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 this 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.197) then
 453                map%source = 'NCEP CDAS CFSV2'
 454              elseif (iprocess.eq.25) then
 455                map%source = 'NCEP SNOW COVER ANALYSIS'
 456              else
 457                map%source = 'unknown model from NCEP'
 458                write (6,*) 'unknown NCEP model, iprocess = ',iprocess
 459              end if
 460            else if (icenter .eq. 57) then
 461              if (iprocess .eq. 87) then
 462                map%source = 'AFWA AGRMET'
 463              else
 464                map%source = 'AFWA'
 465              endif
 466            else if (icenter .eq. 58) then
 467                map%source = 'US Navy FNOC'
 468            else if (icenter .eq. 98) then
 469              map%source = 'ECMWF'
 470            else if (icenter .eq. 34) then
 471              map%source = 'JMA'
 472            else if (icenter .eq. 74 .or. icenter .eq. 75 ) then
 473              map%source = 'UKMO'
 474            else
 475              map%source = 'unknown model and orig center'
 476            end if
 477
 478            !--
 479
 480            ! Store information about the grid on which the data is.
 481            ! This stuff gets stored in the MAP variable, as defined in
 482            ! module GRIDINFO.
 483
 484            map%startloc = 'SWCORNER'
 485
 486            if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
 487               map%igrid = 0
 488               map%nx = gfld%igdtmpl(8)
 489               map%ny = gfld%igdtmpl(9)
 490               map%dx = gfld%igdtmpl(17)
 491               map%dy = gfld%igdtmpl(18)
 492               map%lat1 = gfld%igdtmpl(12)
 493               map%lon1 = gfld%igdtmpl(13)
 494
 495               if ((gfld%igdtmpl(10) .eq. 0).OR.   &
 496                   (gfld%igdtmpl(10) .eq. 255)) THEN
 497           ! Scale lat/lon values to 0-180, default range is 1e6.
 498                 map%lat1 = map%lat1/scale_factor
 499                 map%lon1 = map%lon1/scale_factor
 500           ! Scale dx/dy values to degrees, default range is 1e6.
 501                 map%dx = map%dx/scale_factor
 502                 map%dy = map%dy/scale_factor
 503               else
 504           ! Basic angle and subdivisions are non-zero (not tested)
 505                 map%lat1 = map%lat1 *  &
 506                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 507                 map%lon1 = map%lon1 *  &
 508                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 509                 map%dx = map%dx *  &
 510                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 511                 map%dy = map%dy *  &
 512                            (gfld%igdtmpl(11)/gfld%igdtmpl(10))
 513               endif
 514
 515            elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
 516               map%igrid = 3
 517               map%nx = gfld%igdtmpl(8)
 518               map%ny = gfld%igdtmpl(9)
 519               map%lov = gfld%igdtmpl(14)
 520               map%truelat1 = gfld%igdtmpl(19)
 521               map%truelat2 = gfld%igdtmpl(20)
 522               map%dx = gfld%igdtmpl(15)
 523               map%dy = gfld%igdtmpl(16)
 524               map%lat1 = gfld%igdtmpl(10)
 525               map%lon1 = gfld%igdtmpl(11)
 526
 527            elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
 528               map%igrid = 0
 529               map%nx = gfld%igdtmpl(8)
 530               map%ny = gfld%igdtmpl(9)
 531               map%dx = gfld%igdtmpl(17)
 532               map%dy = gfld%igdtmpl(18) ! ?not in Grid Definition Template 3.40 doc
 533               map%lat1 = gfld%igdtmpl(12)
 534               map%lon1 = gfld%igdtmpl(13)
 535
 536               ! Scale dx/dy values to degrees, default range is 1e6.
 537               if (map%dx.gt.10000) then
 538                  map%dx = map%dx/scale_factor
 539               endif
 540               if (map%dy.gt.10000) then
 541                  map%dy = (map%dy/scale_factor)*(-1)
 542               endif
 543
 544               ! Scale lat/lon values to 0-180, default range is 1e6.
 545               if (map%lat1.ge.scale_factor) then
 546                  map%lat1 = map%lat1/scale_factor
 547               endif
 548               if (map%lon1.ge.scale_factor) then
 549                  map%lon1 = map%lon1/scale_factor
 550               endif
 551            print *,'Gaussian Grid: Dx,Dy,lat,lon',map%dx,map%dy, &
 552              map%lat1,map%lon1
 553
 554            elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
 555               map%igrid = 5
 556               map%nx = gfld%igdtmpl(8)
 557               map%ny = gfld%igdtmpl(9)
 558               !map%lov = gfld%igdtmpl(x) ! ?not in Grid Definition Template 3.20 doc
 559               map%truelat1 = 60.
 560               map%truelat2 = 91.
 561               !map%dx = gfld%igdtmpl(x)
 562               !map%dy = gfld%igdtmpl(x)
 563               map%lat1 = gfld%igdtmpl(10)
 564               map%lon1 = gfld%igdtmpl(11)
 565
 566            else
 567               print*, 'GRIB2 Unknown Projection: ',gfld%igdtnum
 568               print*, 'see Code Table 3.1: Grid Definition Template No'
 569            endif
 570
 571             call gf_free(gfld)
 572          endif
 573
 574          ! ----
 575
 576          ! Continue to unpack GRIB2 field.
 577          if (debug_level .gt. 100) write(6,*) 'numfields = ',numfields
 578          NUM_FIELDS: do n = 1, numfields
 579 ! e.g. U and V would =2, otherwise its usually =1
 580            call gf_getfld(cgrib,lengrib,n,unpack,expand,gfld,ierr)
 581            if (ierr.ne.0) then
 582              write(*,*) ' ERROR extracting field gf_getfld = ',ierr
 583              cycle
 584            endif
 585
 586 ! ------------------------------------
 587          ! Additional print information for developer.
 588          pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),  &
 589                                   gfld%ipdtmpl(2))
 590          if (debug_level .gt. 50 ) then
 591            print *
 592 !          print *,'G2 FIELD ',n
 593            if (n==1) then
 594             write(*,'(/,"GRIB2 SECTION 0 - INDICATOR SECTION:")')
 595             write(*,'(5x,"Discipline"'//pstring) gfld%discipline
 596             write(*,'(5x,"GRIB Edition Number"'//pstring) gfld%version
 597             write(*,'(5x,"GRIB length"'//pstring) lengrib
 598             write(*,'(/,"GRIB2 SECTION 1 - IDENTIFICATION SECTION:")')
 599             write(*,'(5x,"Length of Section"'//pstring) gfld%idsectlen
 600             write(*,'(5x,"Originating Center ID"'//pstring) &
 601                gfld%idsect(1)
 602             write(*,'(5x,"Subcenter ID"'//pstring) gfld%idsect(2)
 603             write(*,'(5x,"GRIB Master Table Version"'//pstring) &
 604             gfld%idsect(3)
 605             write(*,'(5x,"GRIB Local Table Version"'//pstring) &
 606             gfld%idsect(4)
 607             write(*,'(5x,"Significance of Reference Time"'//pstring) &
 608             gfld%idsect(5)
 609             write(*,'(5x,"Year"'//pstring)  gfld%idsect(6)
 610             write(*,'(5x,"Month"'//pstring)  gfld%idsect(7)
 611             write(*,'(5x,"Day"'//pstring)  gfld%idsect(8)
 612             write(*,'(5x,"Hour"'//pstring)  gfld%idsect(9)
 613             write(*,'(5x,"Minute"'//pstring)  gfld%idsect(10)
 614             write(*,'(5x,"Second"'//pstring)  gfld%idsect(11)
 615             write(*,'(5x,"Production Status of data"'//pstring) &
 616               gfld%idsect(12)
 617             write(*,'(5x,"Type of processed data"'//pstring) &
 618               gfld%idsect(13)
 619 !           print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
 620            endif
 621            write(*,'(/,"GRIB2 SECTION 2 - LOCAL SECTION:")')
 622            write(*,'(5x,"Length of Section 2"'//pstring) gfld%locallen
 623            if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
 624            do j = 1, gfld%locallen
 625              write(*,'(5x,"Local value "'//astring) gfld%local(j)
 626            enddo
 627 !             print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
 628            endif
 629            write(*,'(/,"GRIB2 SECTION 3 - GRID DEFINITION SECTION:")')
 630 !          write(*,'(5x,"Length of Section 3"'//pstring) gfld%unknown
 631            write(*,'(5x,"Source of grid definition"'&
 632            //pstring) gfld%griddef
 633            write(*,'(5x,"Number of grid points"'//pstring) gfld%ngrdpts
 634            write(*,'(5x,"Number of octets for addnl points"'//pstring) &
 635            gfld%numoct_opt
 636            write(*,'(5x,"Interpretation list"'//pstring) &
 637            gfld%interp_opt
 638            write(*,'(5x,"Grid Definition Template Number"'//pstring) &
 639            gfld%igdtnum
 640            if (gfld%igdtnum .eq. 0 .or. gfld%igdtnum .eq. 1 .or.  &
 641                gfld%igdtnum .eq. 2 .or. gfld%igdtnum .eq. 3 ) then
 642              if (gfld%igdtnum .eq. 0 ) then
 643                write(*,'(5x,"Lat/Lon or Cylindrical Equidistant Grid")')
 644              else if (gfld%igdtnum .eq. 1 ) then
 645                write(*,'(5x,"Rotated Lat/Lon or Cylind. Equi. Grid")')
 646              else if (gfld%igdtnum .eq. 2 ) then
 647                write(*,'(5x,"Stretched Lat/Lon or Cylind. Equi. Grid")')
 648              else if (gfld%igdtnum .eq. 3 ) then
 649                write(*,'(5x,"Stretched and Rotated Lat/Lon Grid")')
 650              endif
 651              write(*,'(10x,"Shape of the Earth"'//pstring) &
 652                 gfld%igdtmpl(1)
 653              write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
 654                 gfld%igdtmpl(2)
 655              write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
 656                 gfld%igdtmpl(3)
 657              write(*,'(10x,"Scale factor of major axis"'//pstring) &
 658                 gfld%igdtmpl(4)
 659              write(*,'(10x,"Scaled value of major axis"'//pstring) &
 660                 gfld%igdtmpl(5)
 661              write(*,'(10x,"Scale factor of minor axis"'//pstring) &
 662                 gfld%igdtmpl(6)
 663              write(*,'(10x,"Scaled value of minor axis"'//pstring) &
 664                 gfld%igdtmpl(7)
 665              write(*,'(10x,"Ni - points along a parallel"'//pstring) &
 666                  gfld%igdtmpl(8)
 667              write(*,'(10x,"Nj - points along a meridian"'//pstring) &
 668                  gfld%igdtmpl(9)
 669              write(*,'(10x,"Basic angle of initial domain"'//pstring)&
 670                  gfld%igdtmpl(10)
 671              write(*,'(10x,"Subdivisions of basic angle"'//pstring) &
 672                  gfld%igdtmpl(11)
 673              write(*,'(10x,"La1"'//pstring)  gfld%igdtmpl(12)
 674              write(*,'(10x,"Lo1"'//pstring)  gfld%igdtmpl(13)
 675              write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
 676               gfld%igdtmpl(14)
 677              write(*,'(10x,"La2"'//pstring) gfld%igdtmpl(15)
 678              write(*,'(10x,"Lo2"'//pstring) gfld%igdtmpl(16)
 679              write(*,'(10x,"Di - i-dir increment"'//pstring) &
 680                gfld%igdtmpl(17)
 681              write(*,'(10x,"Dj - j-dir increment"'//pstring) &
 682                gfld%igdtmpl(18)
 683              write(*,'(10x,"Scanning mode"'//pstring) &
 684                gfld%igdtmpl(19)
 685              if (gfld%igdtnum .eq. 1) then
 686              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 687                gfld%igdtmpl(20)
 688              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 689                gfld%igdtmpl(21)
 690              write(*,'(10x,"Angle of rotation of projection"'//pstring)&
 691                gfld%igdtmpl(22)
 692              else if (gfld%igdtnum .eq. 2) then
 693              write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
 694                gfld%igdtmpl(20)
 695              write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
 696                gfld%igdtmpl(21)
 697              write(*,'(10x,"Stretching factor"'//pstring) &
 698                gfld%igdtmpl(22)
 699              else if (gfld%igdtnum .eq. 3) then
 700              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 701                gfld%igdtmpl(20)
 702              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 703                gfld%igdtmpl(21)
 704              write(*,'(10x,"Angle of rotation of projection"'//pstring)&
 705                gfld%igdtmpl(22)
 706              write(*,'(10x,"Lat of the pole of stretching "'//pstring)&
 707                gfld%igdtmpl(23)
 708              write(*,'(10x,"Lon of the pole of stretching "'//pstring)&
 709                gfld%igdtmpl(24)
 710              write(*,'(10x,"Stretching factor"'//pstring) &
 711                gfld%igdtmpl(25)
 712              endif
 713            else if (gfld%igdtnum .eq. 10) then
 714              write(*,'(5x,"Mercator Grid")')
 715            else if (gfld%igdtnum .eq. 20 .or. gfld%igdtnum .eq. 30) then
 716              if (gfld%igdtnum .eq. 20) then
 717                write(*,'(5x,"Polar Stereographic Grid")')
 718              else if (gfld%igdtnum .eq. 30) then
 719                write(*,'(5x,"Lambert Conformal Grid")')
 720              endif
 721              write(*,'(10x,"Shape of the Earth"'//pstring) &
 722                 gfld%igdtmpl(1)
 723              write(*,'(10x,"Scale factor of spher. Earth"'//pstring) &
 724                 gfld%igdtmpl(2)
 725              write(*,'(10x,"Scaled value of spher. Earth"'//pstring) &
 726                 gfld%igdtmpl(3)
 727              write(*,'(10x,"Scale factor of major axis"'//pstring) &
 728                 gfld%igdtmpl(4)
 729              write(*,'(10x,"Scaled value of major axis"'//pstring) &
 730                 gfld%igdtmpl(5)
 731              write(*,'(10x,"Scale factor of minor axis"'//pstring) &
 732                 gfld%igdtmpl(6)
 733              write(*,'(10x,"Scaled value of minor axis"'//pstring) &
 734                 gfld%igdtmpl(7)
 735              write(*,'(10x,"Nx"'//pstring) gfld%igdtmpl(8)
 736              write(*,'(10x,"Ny"'//pstring) gfld%igdtmpl(9)
 737              write(*,'(10x,"La1"'//pstring) gfld%igdtmpl(10)
 738              write(*,'(10x,"Lo1"'//pstring) gfld%igdtmpl(11)
 739              write(*,'(10x,"Resolution and Component",t50,":",B14.8)')&
 740               gfld%igdtmpl(12)
 741              write(*,'(10x,"LaD"'//pstring) gfld%igdtmpl(13)
 742              write(*,'(10x,"LoV"'//pstring) gfld%igdtmpl(14)
 743              write(*,'(10x,"Dx"'//pstring) gfld%igdtmpl(15)
 744              write(*,'(10x,"Dy"'//pstring) gfld%igdtmpl(16)
 745              write(*,'(10x,"Projection Center Flag"'//pstring) &
 746                gfld%igdtmpl(17)
 747              write(*,'(10x,"Scanning mode"'//pstring) &
 748                gfld%igdtmpl(18)
 749              if (gfld%igdtnum .eq. 30) then
 750              write(*,'(10x,"Latin 1 "'//pstring) &
 751                gfld%igdtmpl(19)
 752              write(*,'(10x,"Latin 2 "'//pstring) &
 753                gfld%igdtmpl(20)
 754              write(*,'(10x,"Lat of southern pole of project"'//pstring)&
 755                gfld%igdtmpl(21)
 756              write(*,'(10x,"Lon of southern pole of project"'//pstring)&
 757                gfld%igdtmpl(22)
 758              endif
 759            else if (gfld%igdtnum .eq. 40 .or. gfld%igdtnum .eq. 41) then
 760              if (gfld%igdtnum .eq. 40) then
 761                write(*,'(5x,"Gaussian Lat/Lon Grid")')
 762              else if (gfld%igdtnum .eq. 41) then
 763                write(*,'(5x,"Rotated Gaussian Lat/Lon Grid")')
 764              else if (gfld%igdtnum .eq. 42) then
 765                write(*,'(5x,"Stretched Gaussian Lat/Lon Grid")')
 766              else if (gfld%igdtnum .eq. 43) then
 767                write(*,'(5x,"Stretched and Rotated Gaussian Lat/Lon ")')
 768              endif
 769            else
 770            do j = 1, gfld%igdtlen
 771              write(*,'(5x,"Grid Definition Template entry "'//pstring) &
 772              gfld%igdtmpl(j)
 773            enddo
 774            endif
 775 !          print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts, &
 776 !                                 gfld%numoct_opt,gfld%interp_opt,  &
 777 !                                 gfld%igdtnum
 778 !          print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',   &
 779 !                 (gfld%igdtmpl(j),j=1,gfld%igdtlen)
 780            if ( gfld%num_opt .eq. 0 ) then
 781 !            print *,'G2 NO Section 3 List Defining No. of Data Points.'
 782            else
 783              print *,'G2 Section 3 Optional List: ',     &
 784                       (gfld%list_opt(j),j=1,gfld%num_opt)
 785            endif
 786           write(*,'(/,"GRIB2 SECTION 4 - PRODUCT DEFINITION SECTION:")')
 787 !          write(*,'(5x,"Length of Section 4"'//pstring) gfld%unknown
 788            write(*,'(5x,"Product Definition Template Number"'//pstring)&
 789             gfld%ipdtnum
 790             do j = 1, gfld%ipdtlen
 791               write(tmp4,'(i4)') j
 792               string = '(5x,"Template Entry  '//tmp4 // '"'
 793               write(*,string//pstring) gfld%ipdtmpl(j)
 794             enddo
 795 !          print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',   &
 796 !               (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
 797
 798            !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
 799            !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
 800            write(*,'(5x,"Product Abbreviated Name",t50,":",a14)')&
 801            pabbrev
 802 !           print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
 803
 804            if ( gfld%num_coord .eq. 0 ) then
 805 !            print *,'G2 NO Optional Vertical Coordinate List.'
 806            else
 807              print *,'G2 Section 4 Optional Coordinates: ',   &
 808                    (gfld%coord_list(j),j=1,gfld%num_coord)
 809            endif
 810 !          if ( gfld%ibmap .ne. 255 ) then
 811 !             print *,'G2 Num. of Data Points = ',gfld%ndpts,   &
 812 !                  '    with BIT-MAP ',gfld%ibmap
 813 !          else
 814 !             print *,'G2 Num. of Data Points = ',gfld%ndpts,  &
 815 !                     '    NO BIT-MAP '
 816 !          endif
 817          write(*,'(/,"GRIB2 SECTION 5 - DATA REPRESENTATION SECTION:")')
 818          write(*,'(5x,"Data Representation Template Number"'//pstring)&
 819             gfld%idrtnum
 820             do j = 1, gfld%idrtlen
 821               write(tmp4,'(i4)') j
 822               string = '(5x,"Template Entry  '//tmp4 // '"'
 823               write(*,string//pstring) gfld%idrtmpl(j)
 824             enddo
 825 !          print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',  &
 826 !               (gfld%idrtmpl(j),j=1,gfld%idrtlen)
 827
 828 !      if ( gfld%ipdtnum .eq. 0 ) then
 829 !        if (gfld%ipdtmpl(1) .eq. 0 ) then
 830 !          write(6,*) 'Temperature'
 831 !        else if (gfld%ipdtmpl(1) .eq. 1 ) then
 832 !          write(6,*) 'Moisture'
 833 !        else if (gfld%ipdtmpl(1) .eq. 2 ) then
 834 !          write(6,*) 'Momentum'
 835 !        else if (gfld%ipdtmpl(1) .eq. 3 ) then
 836 !         write(6,*) 'Mass'
 837 !       endif
 838 !      endif
 839
 840          write(*,'(/,"GRIB2 SECTION 6 - BIT-MAP SECTION:")')
 841          write(*,'(5x,"Bit-map indicator"'//pstring) &
 842            gfld%ibmap
 843
 844            fldmax=gfld%fld(1)
 845            fldmin=gfld%fld(1)
 846            sum=gfld%fld(1)
 847            do j=2,gfld%ndpts
 848              if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
 849              if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
 850              sum=sum+gfld%fld(j)
 851            enddo ! gfld%ndpts
 852
 853          write(*,'(/,"GRIB2 SECTION 7 - DATA SECTION:")')
 854
 855          write(*,'(5x,"Minimum Data Value "'//rstring)&
 856             fldmin
 857          write(*,'(5x,"Maximum Data Value "'//rstring)&
 858             fldmax
 859 !          print *,'G2 Data Values:'
 860 !          write(*,fmt='("G2 MIN=",f21.8," AVE=",f21.8,    &
 861 !               " MAX=",f21.8)') fldmin,sum/gfld%ndpts,fldmax
 862          if (debug_level .gt. 100 ) then
 863             print*, 'gfld%fld = ',gfld%fld
 864 !           do j=1,gfld%ndpts
 865 !              write(*,*) j, gfld%fld(j)
 866 !           enddo
 867          endif
 868          endif ! Additional Print information
 869 ! ------------------------------------
 870          if ( debug_level .le. 50) then
 871            if(gfld%ipdtmpl(10).eq.100) then    ! pressure level
 872              level=gfld%ipdtmpl(12) *  &
 873                      (10. ** (-1. * gfld%ipdtmpl(11)))
 874            else if(gfld%ipdtmpl(10).eq.101 .or.&    ! sea level, sfc, or trop
 875                 gfld%ipdtmpl(10).eq.1 .or. gfld%ipdtmpl(10).eq.7) then
 876              level = 0
 877            else if(gfld%ipdtmpl(10).eq.106) then  ! below ground sfc is in cm in Vtable
 878              level= 100. * gfld%ipdtmpl(12)*(10.**(-1.*gfld%ipdtmpl(11)))
 879            else
 880              level=gfld%ipdtmpl(12) * 10.** (-1.*gfld%ipdtmpl(11))
 881            endif
 882            if (gfld%ipdtmpl(13) .eq. 255) then
 883              lvl2 = 0
 884            else if(gfld%ipdtmpl(10).eq.106) then   ! below ground sfc is in cm in Vtable
 885              lvl2 = 100. * gfld%ipdtmpl(15) * 10.** (-1.*gfld%ipdtmpl(14))
 886            else
 887              lvl2 = gfld%ipdtmpl(15) * 10.** (-1.*gfld%ipdtmpl(14))
 888            endif
 889 !           Account for multiple forecast hours in one file
 890            if (gfld%ipdtnum.eq.0 .or. gfld%ipdtnum.eq.1 .or. gfld%ipdtnum.eq. 8) then
 891               ! Product Definition Template 4.0, 4.1, 4.8
 892               ! Extract forecast time.
 893              if ( gfld%ipdtmpl(8) .eq. 1 ) then   ! time units are hours
 894                fcst = gfld%ipdtmpl(9)
 895              else if ( gfld%ipdtmpl(8) .eq. 0 ) then  ! minutes
 896                fcst = gfld%ipdtmpl(9) / 60.
 897              else if ( gfld%ipdtmpl(8) .eq. 2 ) then  ! days
 898                fcst = gfld%ipdtmpl(9) * 24.
 899              else
 900                fcst = 999
 901              endif
 902            endif
 903
 904            ! Non-standard Product Definition Templates need to be reported
 905            string = '                      '
 906            if ( gfld%ipdtnum .eq. 8 ) then
 907              string = '  PDT4.8'
 908            else if ( gfld%ipdtnum .eq. 1 ) then
 909              string = '  PDT4.1'
 910            endif
 911            write(6,987) itot,gfld%discipline,gfld%ipdtmpl(1),  &
 912              gfld%ipdtmpl(2),gfld%ipdtmpl(10),int(level),&
 913              lvl2,gfld%ipdtnum,pabbrev,hdate,fcst,string
 914   987     format(2i4,i5,i4,i8,i8,i8,i8,3x,a10,a20,i5.2,a10)
 915
 916          endif
 917
 918          ! Deallocate arrays decoding GRIB2 record.
 919          call gf_free(gfld)
 920
 921          enddo NUM_FIELDS
 922
 923       enddo VERSION ! skgb
 924
 925        if (debug_level .gt. 50) &
 926           print *, 'G2 total number of fields found = ',itot
 927
 928         CALL BACLOSE(junit,IOS)
 929
 930         ireaderr=1
 931       else
 932         print *,'open status failed because',ios
 933         hdate = '9999-99-99_99:99:99'
 934         ireaderr=2
 935       endif ! ireaderr check
 936
 937       END subroutine r_grib2
 938
 939 !*****************************************************************************!
 940 ! Subroutine edition_num                                                         !
 941 !                                                                             !
 942 ! Purpose:                                                                    !
 943 !    Read one record from the input GRIB2 file.  Based on the information in  !
 944 !    the GRIB2 header and the user-defined Vtable, decide whether the field in!
 945 !    the GRIB2 record is one to process or to skip.  If the field is one we   !
 946 !    want to keep, extract the data from the GRIB2 record, and pass the data  !
 947 !    back to the calling routine.                                             !
 948 !                                                                             !
 949 ! Argument list:                                                              !
 950 !    Input:                                                                   !
 951 !       JUNIT   : "Unit Number" to open and read from.  Not really a Fortran  !
 952 !                 unit number, since we do not do Fortran I/O for the GRIB2   !
 953 !                 files.  Nor is it a UNIX File Descriptor returned from a C  !
 954 !                 OPEN statement.  It is really just an array index to the    !
 955 !                 array (IUARR) where the UNIX File Descriptor values are     !
 956 !                 stored.                                                     !
 957 !       GRIB2FILE: File name to open, if it is not already open.              !
 958 !                                                                             !
 959 !    Output:                                                                  !
 960 !       GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2                !
 961 !       IERR     : Error flag: 0 - no error on read from GRIB2 file.          !
 962 !                              1 - Hit the end of the GRIB2 file.             !
 963 !                              2 - The file GRIBFLNM we tried to open does    !
 964 !                                  not exist.                                 !
 965 ! Author: Paula McCaslin                                                      !
 966 ! NOAA/FSL                                                                    !
 967 ! Sept 2004                                                                   !
 968 !*****************************************************************************!
 969
 970       SUBROUTINE edition_num(junit, gribflnm, grib_edition, ireaderr)
 971
 972       use grib_mod
 973       use params
 974
 975       parameter(msk1=32000,msk2=4000)
 976       character(len=1),allocatable,dimension(:) :: cgrib
 977       integer :: listsec0(3)
 978       integer :: listsec1(13)
 979       character(len=*)  :: gribflnm
 980       integer :: lskip, lgrib
 981       integer :: junit
 982       integer :: grib_edition
 983       integer :: i, j, ireaderr
 984       integer :: currlen
 985
 986       character(len=4) :: ctemp
 987       character(len=4),parameter :: grib='GRIB',c7777='7777'
 988
 989 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 990 !  SET ARGUMENTS
 991
 992       itot=0
 993       icount=0
 994       iseek=0
 995       lskip=0
 996       lgrib=0
 997       currlen=0
 998
 999 !/* IOS Return Codes from BACIO:  */
1000 !/*  0    All was well                                   */
1001 !/* -1    Tried to open read only _and_ write only       */
1002 !/* -2    Tried to read and write in the same call       */
1003 !/* -3    Internal failure in name processing            */
1004 !/* -4    Failure in opening file                        */
1005 !/* -5    Tried to read on a write-only file             */
1006 !/* -6    Failed in read to find the 'start' location    */
1007 !/* -7    Tried to write to a read only file             */
1008 !/* -8    Failed in write to find the 'start' location   */
1009 !/* -9    Error in close                                 */
1010 !/* -10   Read or wrote fewer data than requested        */
1011
1012 !if ireaderr =1 we have hit the end of a file.
1013 !if ireaderr =2 we have hit the end of all the files.
1014 !if ireaderr =3 beginning characters 'GRIB' not found
1015
1016 !     write(6,*) 'junit = ',junit,' gribflnm = ',gribflnm
1017
1018       ! Open a byte-addressable file.
1019       CALL BAOPENR(junit,gribflnm,IOS)      ! from w3lib
1020 !     write(6,*) 'ios = ',ios
1021       if (ios.eq.0) then
1022
1023          ! Search opend file for the next GRIB2 messege (record).
1024          call skgb(junit,iseek,msk1,lskip,lgrib)
1025
1026          ! Check for EOF, or problem
1027          if (lgrib.eq.0) then
1028             write(*,'("\n\tThere is a problem with the input file.")')
1029             write(*,'("\tPerhaps it is not a Grib2 file?\n")')
1030             STOP "Grib2 file or date problem, stopping in edition_num."
1031          endif
1032
1033          ! Check size, if needed allocate more memory.
1034          if (lgrib.gt.currlen) then
1035             if (allocated(cgrib)) deallocate(cgrib)
1036             allocate(cgrib(lgrib),stat=is)
1037             currlen=lgrib
1038          endif
1039
1040          ! Read a given number of bytes from unblocked file.
1041          call baread(junit,lskip,lgrib,lengrib,cgrib)
1042
1043          ! Check for beginning of GRIB message in the first 100 bytes
1044          istart=0
1045          do j=1,100
1046             ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
1047             if (ctemp.eq.grib ) then
1048               istart=j
1049               exit
1050             endif
1051          enddo
1052          if (istart.eq.0) then
1053             ireaderr=3
1054             print*, "The beginning 4 characters >GRIB< were not found."
1055          endif
1056
1057          ! Unpack Section 0 - Indicator Section to extract GRIB edition field
1058          iofst=8*(istart+5)
1059          call gbyte(cgrib,discipline,iofst,8)     ! Discipline
1060          iofst=iofst+8
1061          call gbyte(cgrib,grib_edition,iofst,8)   ! GRIB edition number
1062
1063          print *, 'ungrib - grib edition num',  grib_edition
1064          CALL BACLOSE(junit,IOS)
1065          ireaderr=1
1066       else if (ios .eq. -4) then
1067         print *,'edition_num: unable to open ',gribflnm
1068         stop 'edition_num'
1069       else
1070          print *,'edition_num: open status failed because',ios,gribflnm
1071          ireaderr=2
1072       endif ! ireaderr check
1073
1074       END subroutine edition_num
1075 subroutine parse_args(err, a1, h1, i1, l1, a2, h2, i2, l2, a3, h3, i3, l3, &
1076      hlast)
1077   integer :: err
1078   character(len=*) , optional :: a1, a2, a3
1079   character(len=*), optional :: h1, h2, h3
1080   integer , optional :: i1, i2, i3
1081   logical, optional :: l1, l2, l3
1082   character(len=*), optional :: hlast
1083
1084   character(len=100) :: hold
1085   integer :: ioff = 0
1086
1087   if (present(hlast)) then
1088      ioff = -1
1089   endif
1090
1091   err = 0
1092
1093   narg = iargc()
1094   numarg = narg + ioff
1095
1096   i = 1
1097   LOOP : do while ( i <= numarg)
1098
1099      ierr = 1
1100      if (present(i1)) then
1101         call checkiarg(i, a1, i1, ierr)
1102      elseif (present(h1)) then
1103         call checkharg(i, a1, h1, ierr)
1104      elseif (present(l1)) then
1105         call checklarg(i, a1, l1, ierr)
1106      endif
1107      if (ierr.eq.0) cycle LOOP
1108
1109      if (present(i2)) then
1110         call checkiarg(i, a2, i2, ierr)
1111      elseif (present(h2)) then
1112         call checkharg(i, a2, h2, ierr)
1113      elseif (present(l2)) then
1114         call checklarg(i, a2, l2, ierr)
1115      endif
1116      if (ierr.eq.0) cycle LOOP
1117
1118      if (present(i3)) then
1119         call checkiarg(i, a3, i3, ierr)
1120      elseif (present(h3)) then
1121         call checkharg(i, a3, h3, ierr)
1122      elseif (present(l3)) then
1123         call checklarg(i, a3, l3, ierr)
1124      endif
1125      if (ierr.eq.0) cycle LOOP
1126
1127      err = 1
1128      call getarg(1, hold)
1129      write(*, '("arg = ", A)') trim(hold)
1130
1131      exit LOOP
1132
1133   enddo LOOP
1134
1135   if (present(hlast)) then
1136      if (narg.eq.0) then
1137         err = 1
1138      else
1139         call getarg(narg, hlast)
1140      endif
1141   endif
1142
1143 contains
1144   subroutine checkiarg(c, a, i, ierr)
1145     integer :: c
1146     character(len=*) :: a
1147     integer :: i
1148
1149     character(len=100) :: hold
1150     ierr = 1
1151
1152     call getarg(c, hold)
1153
1154     if ('-'//a.eq.trim(hold)) then
1155        c = c + 1
1156        call getarg(c, hold)
1157        read(hold, *) i
1158        c = c + 1
1159        ierr = 0
1160     elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1161        hold = hold(len_trim(a)+2: len(hold))
1162        read(hold, *) i
1163        c = c + 1
1164        ierr = 0
1165     endif
1166
1167   end subroutine checkiarg
1168   subroutine checkharg(c, a, h, ierr)
1169     integer :: c
1170     character(len=*) :: a
1171     character(len=*) :: h
1172
1173     character(len=100) :: hold
1174     ierr = 1
1175
1176     call getarg(c, hold)
1177
1178     if ('-'//a.eq.trim(hold)) then
1179        c = c + 1
1180        call getarg(c, hold)
1181        h = trim(hold)
1182        c = c + 1
1183        ierr = 0
1184     elseif ('-'//a .eq. hold(1:len_trim(a)+1)) then
1185        hold = hold(len_trim(a)+2: len(hold))
1186        h = trim(hold)
1187        c = c + 1
1188        ierr = 0
1189     endif
1190
1191   end subroutine checkharg
1192
1193   subroutine checklarg(c, a, l, ierr)
1194     integer :: c
1195     character(len=*) :: a
1196     logical :: l
1197
1198     character(len=100) :: hold
1199     ierr = 1
1200
1201     call getarg(c, hold)
1202     if ('-'//a.eq.trim(hold)) then
1203        l = .TRUE.
1204        c = c + 1
1205        ierr = 0
1206     endif
1207
1208   end subroutine checklarg
1209
1210 end subroutine parse_args
1211