ungrib/src/g2print.F90

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