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