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