WPS/ungrib/src/rd_grib2.F90

   1 !*****************************************************************************!
   2 ! Subroutine RD_GRIB2                                                         !
   3 !                                                                             !
   4 ! Purpose:                                                                    !
   5 !    Read one record from the input GRIB2 file.  Based on the information in  !
   6 !    the GRIB2 header and the user-defined Vtable, decide whether the field in!
   7 !    the GRIB2 record is one to process or to skip.  If the field is one we   !
   8 !    want to keep, extract the data from the GRIB2 record, and store the data !
   9 !    in the ungrib memory structure.                                          !
  10 !                                                                             !
  11 ! Argument list:                                                              !
  12 !    Input:                                                                   !
  13 !       junit   : "Unit Number" to open and read from.  Not really a Fortran  !
  14 !                 unit number, since we do not do Fortran I/O for the GRIB2   !
  15 !                 files.  Nor is it a UNIX File Descriptor returned from a C  !
  16 !                 OPEN statement.  It is really just an array index to the    !
  17 !                 array (IUARR) where the UNIX File Descriptor values are     !
  18 !                 stored.                                                     !
  19 !       gribflnm     : File name to open, if it is not already open.          !
  20 !       debug_level  : Integer for various amounts of printout.               !
  21 !                                                                             !
  22 !    Output:                                                                  !
  23 !                                                                             !
  24 !       hdate        : The (up to)19-character date of the field to process.  !
  25 !       grib_edition : Version of the gribfile (1 or 2)                       !
  26 !       ireaderr     : Error flag: 0 - no error on read from GRIB2 file.      !
  27 !                              1 - Hit the end of the GRIB2 file.             !
  28 !                              2 - The file GRIBFLNM we tried to open does    !
  29 !                                  not exist.                                 !
  30 !                                                                             !
  31 !                                                                             !
  32 ! Author: Paula McCaslin, NOAA/FSL,   Sept 2004                               !
  33 ! Code is based on code developed by Steve Gilbert NCEP & Kevin Manning NCAR  !
  34 ! Adapted for WPS: Jim Bresch, NCAR/MMM. Sept 2006                            !
  35 !*****************************************************************************!
  36
  37       SUBROUTINE rd_grib2(junit, gribflnm, hdate,
  38      &  grib_edition, ireaderr, debug_level)
  39
  40       use grib_mod
  41       use params
  42       use table          ! Included to define g2code
  43       use gridinfo       ! Included to define map%
  44       use storage_module ! Included sub put_storage
  45       use module_debug
  46
  47       real, allocatable, dimension(:) :: hold_array
  48       parameter(msk1=32000,msk2=4000)
  49       character(len=1),allocatable,dimension(:) :: cgrib
  50       integer :: listsec0(3)
  51       integer :: listsec1(13)
  52       integer year, month, day, hour, minute, second, fcst
  53       character(len=*)  :: gribflnm
  54       character(len=*)  :: hdate
  55       character(len=8)  :: pabbrev
  56       character(len=20) :: labbrev
  57       character(len=80) :: tabbrev
  58       integer :: lskip, lgrib
  59       integer :: junit, itot, icount, iseek
  60       integer :: grib_edition
  61       integer :: i, j, ireaderr, ith , debug_level
  62       integer :: currlen
  63       logical :: unpack, expand
  64       type(gribfield) :: gfld
  65       ! For subroutine put_storage
  66       real :: level
  67       real :: scale_factor
  68       integer :: iplvl
  69       character (len=9) :: my_field
  70       character (len=8) :: tmp8
  71       ! For subroutine outout
  72       integer , parameter :: maxlvl = 100
  73       real , dimension(maxlvl) :: plvl
  74       integer :: nlvl
  75       integer , dimension(maxlvl) :: level_array
  76
  77 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  78 C  SET ARGUMENTS
  79
  80       unpack=.true.
  81       expand=.true.
  82       hdate = '0000-00-00_00:00:00'
  83       ierr=0
  84       itot=0
  85       icount=0
  86       iseek=0
  87       lskip=0
  88       lgrib=0
  89       currlen=0
  90       ith=1
  91       scale_factor = 1e6
  92       call mprintf(.true.,DEBUG,"Begin rd_grib2", newline=.true.)
  93
  94 !/* IOS Return Codes from BACIO:  */
  95 !/*  0    All was well                                   */
  96 !/* -1    Tried to open read only _and_ write only       */
  97 !/* -2    Tried to read and write in the same call       */
  98 !/* -3    Internal failure in name processing            */
  99 !/* -4    Failure in opening file                        */
 100 !/* -5    Tried to read on a write-only file             */
 101 !/* -6    Failed in read to find the 'start' location    */
 102 !/* -7    Tried to write to a read only file             */
 103 !/* -8    Failed in write to find the 'start' location   */
 104 !/* -9    Error in close                                 */
 105 !/* -10   Read or wrote fewer data than requested        */
 106
 107 !if ireaderr =1 we have hit the end of a file.
 108 !if ireaderr =2 we have hit the end of all the files.
 109
 110
 111       ! Open a byte-addressable file.
 112       CALL BAOPENR(junit,gribflnm,IOS)
 113       if (ios.eq.0) then
 114       VERSION: do
 115
 116          ! Search opend file for the next GRIB2 messege (record).
 117          call skgb(junit,iseek,msk1,lskip,lgrib)
 118
 119          ! Check for EOF, or problem
 120          if (lgrib.eq.0) then
 121             exit
 122          endif
 123
 124          ! Check size, if needed allocate more memory.
 125          if (lgrib.gt.currlen) then
 126             if (allocated(cgrib)) deallocate(cgrib)
 127             allocate(cgrib(lgrib),stat=is)
 128             !print *,'G2 allocate(cgrib(lgrib)) status: ',IS
 129             currlen=lgrib
 130          endif
 131
 132          ! Read a given number of bytes from unblocked file.
 133          call baread(junit,lskip,lgrib,lengrib,cgrib)
 134
 135          call mprintf ((lgrib.ne.lengrib),ERROR,
 136      &    "rd_grib2: IO Error. %i .ne. %i ", newline=.true.,
 137      &    i1=lgrib,i2=lengrib)
 138
 139          iseek=lskip+lgrib
 140          icount=icount+1
 141
 142          call mprintf (.true.,DEBUG,
 143      &     "G2 GRIB MESSAGE  %i starts at %i ", newline=.true.,
 144      &      i1=icount,i2=lskip+1)
 145
 146          ! Unpack GRIB2 field
 147          call gb_info(cgrib,lengrib,listsec0,listsec1,
 148      &                numfields,numlocal,maxlocal,ierr)
 149          call mprintf((ierr.ne.0),ERROR,
 150      &     " ERROR querying GRIB2 message = %i",newline=.true.,i1=ierr)
 151          itot=itot+numfields
 152
 153          grib_edition=listsec0(2)
 154          if (grib_edition.ne.2) then
 155               exit VERSION
 156          endif
 157
 158          ! Additional print statments for developer.
 159 !MGD         if ( debug_level .GT. 100 ) then
 160 !MGD     print *,'G2 SECTION 0: ',(listsec0(j),j=1,3)
 161 !MGD         print *,'G2 SECTION 1: ',(listsec1(j),j=1,13)
 162 !MGD         print *,'G2 Contains ',numlocal,' Local Sections ',
 163 !MGD     &           ' and ',numfields,' data fields.'
 164 !MGD         endif
 165
 166
 167          ! ----
 168          ! Once per file fill in date, model and projection values.
 169
 170          if (lskip.lt.10) then
 171
 172            ! Build the 19-character date string, based on GRIB2 header date
 173            ! and time information, including forecast time information:
 174
 175            n=1
 176            call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
 177            year  =gfld%idsect(6)     !(FOUR-DIGIT) YEAR OF THE DATA
 178            month =gfld%idsect(7)     ! MONTH OF THE DATA
 179            day   =gfld%idsect(8)     ! DAY OF THE DATA
 180            hour  =gfld%idsect(9)     ! HOUR OF THE DATA
 181            minute=gfld%idsect(10)    ! MINUTE OF THE DATA
 182            second=gfld%idsect(11)    ! SECOND OF THE DATA
 183
 184            fcst = 0
 185            ! Parse the forecast time info from Sect 4.
 186            if (gfld%ipdtnum.eq.0) then  ! Product Definition Template 4.0
 187
 188              ! Extract forecast time.
 189              fcst = gfld%ipdtmpl(9)
 190
 191            endif
 192
 193            ! Compute valid time.
 194
 195            !print *, 'ymd',gfld%idsect(6),gfld%idsect(7),gfld%idsect(8)
 196            !print *, 'hhmm  ',gfld%idsect(9),gfld%idsect(10)
 197
 198            call build_hdate(hdate,year,month,day,hour,minute,second)
 199            call mprintf(.true.,DEBUG,"G2 hdate = %s ", newline=.true.,
 200      &                  s1=hdate)
 201            call geth_newdate(hdate,hdate,3600*fcst)
 202            call mprintf(.true.,DEBUG,"G2 hdate (fcst?) = %s ",
 203      &                  newline=.true., s1=hdate)
 204
 205            !--
 206
 207            ! Indicator of the source (center) of the data.
 208            icenter = gfld%idsect(1)
 209
 210            ! Indicator of model (or whatever) which generated the data.
 211            iprocess = gfld%ipdtmpl(5)
 212
 213
 214            if (icenter.eq.7) then
 215              if (iprocess.eq.83 .or. iprocess.eq.84) then
 216                map%source = 'NCEP NAM Model'
 217              elseif (iprocess.eq.81) then
 218                map%source = 'NCEP GFS Model'
 219              elseif (iprocess.eq.96) then
 220                map%source = 'NCEP GFS Model'
 221              elseif (iprocess.eq.109) then
 222                map%source = 'NCEP RTMA'
 223              elseif (iprocess.eq.105) then
 224                map%source = 'NCEP RUC Model'
 225              elseif (iprocess.eq.140) then
 226                map%source = 'NCEP NARR'
 227              elseif (iprocess.eq.44) then
 228                map%source = 'NCEP SST Analysis'
 229              else
 230                map%source = 'unknown model from NCEP'
 231                call mprintf(.true.,STDOUT,
 232      &            "unknown model from NCEP %i ",newline=.true.,
 233      &            i1=iprocess)
 234                call mprintf(.true.,LOGFILE,
 235      &            "unknown model from NCEP %i ",newline=.true.,
 236      &            i1=iprocess)
 237              end if
 238            else if (icenter .eq. 57) then
 239              if (iprocess .eq. 87) then
 240                map%source = 'AFWA AGRMET'
 241              else
 242                map%source = 'AFWA'
 243              endif
 244            else if (icenter .eq. 98) then
 245              map%source = 'ECMWF'
 246            else
 247              map%source = 'unknown model and orig center'
 248            end if
 249
 250            !--
 251
 252            ! Store information about the grid on which the data is.
 253            ! This stuff gets stored in the MAP variable, as defined in
 254            ! module GRIDINFO.
 255
 256            map%startloc = 'SWCORNER'
 257            map%grid_wind = .true.
 258
 259            if (gfld%igdtnum.eq.0) then ! Lat/Lon grid aka Cylindrical Equidistant
 260               map%igrid = 0
 261               map%nx = gfld%igdtmpl(8)
 262               map%ny = gfld%igdtmpl(9)
 263               map%dx = gfld%igdtmpl(17)
 264               map%dy = gfld%igdtmpl(18)
 265               map%lat1 = gfld%igdtmpl(12)
 266               map%lon1 = gfld%igdtmpl(13)
 267               write(tmp8,'(b8.8)') gfld%igdtmpl(14)
 268               if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 269               map%r_earth = earth_radius (gfld%igdtmpl(1))
 270
 271               ! Scale dx/dy values to degrees, default range is 1e6.
 272               if (map%dx.gt.10000) then
 273                  map%dx = map%dx/scale_factor
 274               endif
 275               if (map%dy.gt.10000) then
 276                  map%dy = map%dy/scale_factor
 277               endif
 278
 279               ! Scale lat/lon values to 0-180, default range is 1e6.
 280               if (map%lat1.ge.scale_factor) then
 281                  map%lat1 = map%lat1/scale_factor
 282               endif
 283               if (map%lon1.ge.scale_factor) then
 284                  map%lon1 = map%lon1/scale_factor
 285               endif
 286
 287 ! The following is needed for NCEP GFS, 0.5 degree output. The j-scan is in the -y direction.
 288 ! In WPS, the sign of dy indicates the direction of the scan.
 289               write(tmp8,'(b8.8)') gfld%igdtmpl(19)
 290               read(tmp8,'(1x,i1)') jscan
 291               if ( jscan .eq. 0 .and. map%dy .gt. 0. ) then
 292                 map%dy = -1. * map%dy
 293               endif
 294 !             if ( map%lat1 .gt. gfld%igdtmpl(15) .and.
 295 !    &               map%dy .gt. 0. ) then
 296 !               map%dy = -1. * map%dy
 297 !               write(6,*) 'Resetting map%dy for iprocess = ',iprocess
 298 !             endif
 299
 300            elseif (gfld%igdtnum.eq.20) then ! Polar-Stereographic Grid.
 301               map%igrid = 5
 302               map%nx = gfld%igdtmpl(8)
 303               map%ny = gfld%igdtmpl(9)
 304               map%lov = gfld%igdtmpl(14) / scale_factor
 305               map%truelat1 = 60.
 306               map%truelat2 = 91.
 307               map%dx = gfld%igdtmpl(15) / scale_factor
 308               map%dy = gfld%igdtmpl(16) / scale_factor
 309               map%lat1 = gfld%igdtmpl(10) / scale_factor
 310               map%lon1 = gfld%igdtmpl(11) / scale_factor
 311               write(tmp8,'(b8.8)') gfld%igdtmpl(12)
 312               if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 313               map%r_earth = earth_radius (gfld%igdtmpl(1))
 314
 315            elseif (gfld%igdtnum.eq.30) then ! Lambert Conformal Grid
 316               map%igrid = 3
 317               map%nx = gfld%igdtmpl(8)
 318               map%ny = gfld%igdtmpl(9)
 319               map%lov = gfld%igdtmpl(14) / scale_factor
 320               map%truelat1 = gfld%igdtmpl(19) / scale_factor
 321               map%truelat2 = gfld%igdtmpl(20) / scale_factor
 322               map%dx = gfld%igdtmpl(15) / scale_factor
 323               map%dy = gfld%igdtmpl(16) / scale_factor
 324               map%lat1 = gfld%igdtmpl(10) / scale_factor
 325               map%lon1 = gfld%igdtmpl(11) / scale_factor
 326               write(tmp8,'(b8.8)') gfld%igdtmpl(12)
 327               if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 328               map%r_earth = earth_radius (gfld%igdtmpl(1))
 329
 330            elseif(gfld%igdtnum.eq.40) then ! Gaussian Grid (we will call it lat/lon)
 331               map%igrid = 4
 332               map%nx = gfld%igdtmpl(8)     ! Ni - # of points along a parallel
 333               map%ny = gfld%igdtmpl(9)     ! Nj - # of points along meridian
 334               map%dx = gfld%igdtmpl(17)    ! Di - i direction increment
 335               map%dy = gfld%igdtmpl(18)    ! N - # of parallels between pole and equator
 336               map%lat1 = gfld%igdtmpl(12)  ! La1 - lat of 1st grid point
 337               map%lon1 = gfld%igdtmpl(13)  ! Lo1 - lon of 1st grid point
 338               write(tmp8,'(b8.8)') gfld%igdtmpl(14)  ! resolution/component flag
 339               if (tmp8(5:5) .eq. '0') map%grid_wind = .false.
 340               map%r_earth = earth_radius (gfld%igdtmpl(1))
 341
 342               ! Scale dx/dy values to degrees, default range is 1e6.
 343               if (map%dx.gt.10000) then
 344                  map%dx = map%dx/scale_factor
 345               endif
 346               if (map%dy.gt.10000) then
 347                  map%dy = (map%dy/scale_factor)*(-1)
 348               endif
 349
 350               ! Scale lat/lon values to 0-180, default range is 1e6.
 351               if (map%lat1.ge.scale_factor) then
 352                  map%lat1 = map%lat1/scale_factor
 353               endif
 354               if (map%lon1.ge.scale_factor) then
 355                  map%lon1 = map%lon1/scale_factor
 356               endif
 357               if ( debug_level .gt. 2 ) then
 358               call mprintf(.true.,DEBUG,
 359      &     "Gaussian Grid: Dx,Dy,lat,lon,nlats %f %f %f %f %i ",
 360      &     newline=.true.,f1=map%dx,f2=map%dy,f3=map%lat1,f4=map%lon1,
 361      &     i1=nint(map%dy))
 362               end if
 363
 364            else
 365               call mprintf(.true.,STDOUT,"GRIB2 Unknown Projection: %i",
 366      &          newline=.true.,i1=gfld%igdtnum)
 367               call mprintf(.true.,STDOUT,
 368      &          "see Code Table 3.1: Grid Definition Template Number",
 369      &          newline=.true.)
 370               call mprintf(.true.,LOGFILE,
 371      &          "GRIB2 Unknown Projection: %i",
 372      &          newline=.true.,i1=gfld%igdtnum)
 373               call mprintf(.true.,LOGFILE,
 374      &          "see Code Table 3.1: Grid Definition Template Number",
 375      &          newline=.true.)
 376            endif
 377
 378            if (icenter.eq.7) then
 379              call ncep_grid_num (gfld%igdtnum)
 380            endif
 381
 382            ! Deallocate arrays decoding GRIB2 record.
 383            call gf_free(gfld)
 384          endif
 385
 386          ! ----
 387
 388          ! Continue to unpack GRIB2 field.
 389          do n=1,numfields ! e.g. U and V would =2, otherwise its usually =1
 390            call gf_getfld(cgrib,lengrib,n,.FALSE.,expand,gfld,ierr)
 391            if (ierr.ne.0) then
 392              write(*,*) ' ERROR extracting field gf_getfld = ',ierr
 393              cycle
 394            endif
 395
 396 ! ------------------------------------
 397          ! Additional print information for developer.
 398          if ( debug_level .GT. 1000 ) then
 399 !MGD           print *
 400 !MGD           print *,'G2 FIELD ',n
 401 !MGD           if (n==1) then
 402 !MGD            print *,'G2 SECTION 0: ',gfld%discipline,gfld%version
 403 !MGD            print *,'G2 SECTION 1: ',(gfld%idsect(j),j=1,gfld%idsectlen)
 404 !MGD           endif
 405 !MGD           if ( associated(gfld%local).AND.gfld%locallen.gt.0 ) then
 406 !MGD              print *,'G2 SECTION 2: ',(gfld%local(j),j=1,gfld%locallen)
 407 !MGD           endif
 408 !MGD           print *,'G2 SECTION 3: ',gfld%griddef,gfld%ngrdpts,
 409 !MGD     &                            gfld%numoct_opt,gfld%interp_opt,
 410 !MGD     &                            gfld%igdtnum
 411 !MGD           print *,'G2 GRID TEMPLATE 3.',gfld%igdtnum,': ',
 412 !MGD     &            (gfld%igdtmpl(j),j=1,gfld%igdtlen)
 413 !MGD           if ( gfld%num_opt .eq. 0 ) then
 414 !MGD             print *,'G2 NO Section 3 List Defining No. of Data Points.'
 415 !MGD           else
 416 !MGD             print *,'G2 Section 3 Optional List: ',
 417 !MGD     &                (gfld%list_opt(j),j=1,gfld%num_opt)
 418 !MGD           endif
 419 !MGD           print *,'G2 PRODUCT TEMPLATE 4.',gfld%ipdtnum,': ',
 420 !MGD     &          (gfld%ipdtmpl(j),j=1,gfld%ipdtlen)
 421
 422            pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
 423      &                              gfld%ipdtmpl(2))
 424            !call prlevel(gfld%ipdtnum,gfld%ipdtmpl,labbrev)
 425            !call prvtime(gfld%ipdtnum,gfld%ipdtmpl,listsec1,tabbrev)
 426 !MGD            print *,'G2 TEXT: ',pabbrev,trim(labbrev)," ",trim(tabbrev)
 427
 428 !MGD           if ( gfld%num_coord .eq. 0 ) then
 429 !MGD             print *,'G2 NO Optional Vertical Coordinate List.'
 430 !MGD           else
 431 !MGD             print *,'G2 Section 4 Optional Coordinates: ',
 432 !MGD     &             (gfld%coord_list(j),j=1,gfld%num_coord)
 433 !MGD           endif
 434            if ( gfld%ibmap .ne. 255 ) then
 435               call mprintf(.true.,DEBUG,
 436      &             'G2 Num. of Data Points = %i with BIT-MAP %i',
 437      &             newline=.true., i1=gfld%ndpts, i2=gfld%ibmap)
 438            else
 439               call mprintf(.true.,DEBUG,
 440      &             'G2 Num. of Data Points = %i NO BIT-MAP',
 441      &             newline=.true., i1=gfld%ndpts)
 442            endif
 443 !MGD           print *,'G2 DRS TEMPLATE 5.',gfld%idrtnum,': ',
 444 !MGD     &          (gfld%idrtmpl(j),j=1,gfld%idrtlen)
 445            fldmax=gfld%fld(1)
 446            fldmin=gfld%fld(1)
 447            sum=gfld%fld(1)
 448            do j=2,gfld%ndpts
 449              if (gfld%fld(j).gt.fldmax) fldmax=gfld%fld(j)
 450              if (gfld%fld(j).lt.fldmin) fldmin=gfld%fld(j)
 451              sum=sum+gfld%fld(j)
 452            enddo ! gfld%ndpts
 453
 454 !MGD           print *,'G2 Data Values:'
 455            call mprintf(.true.,DEBUG,'G2 MIN=%f AVE=%f MAX=%f',
 456      &         newline=.true., f1=fldmin, f2=sum/gfld%ndpts, f3=fldmax)
 457            !do j=1,gfld%ndpts\20
 458            !   write(*,*) j, gfld%fld(j)
 459            !enddo
 460          endif ! Additional Print information
 461 ! ------------------------------------
 462
 463 !         do i = 1, maxvar
 464 !           write(6,'(a10,4i8)') namvar(i),(g2code(j,i),j=1,4)
 465 !         enddo
 466 !MGD      if (debug_level .gt. 50) then
 467 !MGD          write(6,*) 'looking for ',gfld%discipline,gfld%ipdtmpl(1),
 468 !MGD     &       gfld%ipdtmpl(2),gfld%ipdtmpl(10)
 469 !MGD          endif
 470
 471          ! Test this data record against list of desired variables
 472          ! found in Vtable.
 473          ! ----
 474          MATCH_LOOP: do i=1,maxvar ! Max variables found in Vtable,
 475                                    ! maxvar is defined in table.mod
 476
 477            if ( gfld%discipline .eq. g2code(1,i) .and.   !Discipline
 478      &          gfld%ipdtmpl(1) .eq. g2code(2,i) .and.   !Category
 479      &          gfld%ipdtmpl(2) .eq. g2code(3,i) .and.   !Parameter
 480      &          gfld%ipdtmpl(10) .eq. g2code(4,i)) then  !Elevation
 481
 482             call gf_free(gfld)
 483             call gf_getfld(cgrib,lengrib,n,.TRUE.,expand,gfld,ierr)
 484             pabbrev=param_get_abbrev(gfld%discipline,gfld%ipdtmpl(1),
 485      &                               gfld%ipdtmpl(2))
 486
 487               !my_field (e.g. RH, TMP, similar to, but not the same as pabbrev)
 488               my_field=namvar(i)
 489
 490 !MGD    if (debug_level .gt. 50) then
 491 !MGD     write(6,*) 'namvar(i) = ',namvar(i),' pabbrev = ',pabbrev
 492 !MGD     write(6,*) 'Parameter = ',gfld%ipdtmpl(2)
 493 !MGD    endif
 494
 495
 496 ! need to match up soil levels with those requested.
 497 ! For the Vtable levels, -88 = all levels, -99 = missing. The units
 498 ! vary depending on the level code (e.g. 106 = cm, 103 = m).
 499               if ( gfld%ipdtmpl(10) .eq. 106 ) then
 500                 TMP8LOOP: do j = 1, maxvar
 501                   if ((g2code(4,j) .eq. 106) .and.
 502      &               (gfld%ipdtmpl(2) .eq. g2code(3,j)) .and.
 503      &               (gfld%ipdtmpl(12) .eq. level1(j)) .and.
 504      &               ((gfld%ipdtmpl(15) .eq. level2(j)) .or.
 505      &                                   (level2(j) .le. -88))) then
 506                     my_field = namvar(j)
 507                     exit TMP8LOOP
 508                   endif
 509                 enddo TMP8LOOP
 510                 if (j .gt. maxvar ) then
 511                   write(6,'(a,i6,a,i6,a)') 'Subsoil level ',
 512      &               gfld%ipdtmpl(12),' to ',gfld%ipdtmpl(15),
 513      &           ' in the GRIB2 file, was not found in the Vtable'
 514                   cycle MATCH_LOOP
 515                 endif
 516 !MGD         if (debug_level .gt. 50) write(6,*) 'my_field is now ',my_field
 517               endif
 518
 519               ! Level (eg. 10000 mb)
 520               if(gfld%ipdtmpl(10).eq.100) then
 521                  ! Pressure level (range from 1000mb to 0mb)
 522                  level=gfld%ipdtmpl(12)
 523               elseif(gfld%ipdtmpl(10).eq.105) then
 524                  ! Hybrid level (range from 1 to N)
 525                  level=gfld%ipdtmpl(12)
 526               elseif(gfld%ipdtmpl(10).eq.104) then
 527                  ! Sigma level (range from 10000 to 0)
 528                  level=gfld%ipdtmpl(12)
 529               elseif(gfld%ipdtmpl(10).eq.101) then
 530                  ! MSL
 531                  level=201300.
 532               elseif(gfld%ipdtmpl(10).eq.106.or.
 533      &               gfld%ipdtmpl(10).eq.1) then
 534                  ! Misc near ground/surface levels
 535                  level=200100.
 536               else
 537                  ! Misc near ground/surface levels
 538                  level=200100.
 539               endif
 540               iplvl = int(level)
 541
 542               ! Store the unpacked 2D slab from the GRIB2 record
 543               allocate(hold_array(gfld%ngrdpts))
 544               do j=1,gfld%ngrdpts
 545                  hold_array(j)=gfld%fld(j)
 546               enddo
 547
 548 !   Some grids need to be reordered. Until we get an example, this is
 549 !   a placeholder
 550 !             call reorder_it (hold_array, map%nx, map%ny, map%dx,
 551 !    &                 map%dy, iorder)
 552
 553               ! When we have reached this point, we have a data array ARRAY
 554               ! which has some data we want to save, with field name FIELD
 555               ! at pressure level LEVEL (Pa).  Dimensions of this data are
 556               ! map%nx and map%ny.  Put this data into storage.
 557
 558               !print *,'call put_storage',iplvl,my_field,hold_array(55),ith
 559               !e.g. call put_storage(200100, 'RH', my_field, 1, ith)
 560               call put_storage(iplvl,my_field,
 561      &           reshape(hold_array(1:map%nx*map%ny),
 562      &           (/map%nx, map%ny/)), map%nx,map%ny)
 563               deallocate(hold_array)
 564
 565               ! If Specific Humidity is present on hybrid levels AND
 566               ! upper-air RH is missing, see if we can compute RH from
 567               ! Specific Humidity.
 568               if (.not. is_there(iplvl, 'RH') .and.
 569      &            is_there(iplvl, 'SH') .and.
 570      &            is_there(iplvl, 'TT') .and.
 571      &            is_there(iplvl, 'P')) then
 572                   call g2_compute_rh_spechumd_upa(map%nx,map%ny,iplvl)
 573                  !call llstor_remove(iplvl, 'SH') !We are done with SH
 574               endif
 575
 576               ith=ith+1
 577               exit MATCH_LOOP
 578
 579            endif ! Selected param.
 580
 581
 582          enddo MATCH_LOOP
 583
 584          ! Deallocate arrays decoding GRIB2 record.
 585          call gf_free(gfld)
 586
 587          enddo ! 1,numfields
 588
 589
 590       enddo VERSION ! skgb
 591
 592
 593        if ( debug_level .gt. 100 ) then
 594          call mprintf (.true.,DEBUG,
 595      &   "G2 total number of fields found = %i ",newline=.true.,i1=itot)
 596        end if
 597
 598        CALL BACLOSE(junit,IOS)
 599
 600        nullify(gfld%local)            ! must be nullified before opening next file
 601        ireaderr=1
 602       else
 603         call mprintf (.true.,DEBUG,"open status failed because %i ",
 604      &                newline=.true., i1=ios)
 605         hdate = '9999-99-99_99:99:99'
 606         ireaderr=2
 607       endif ! ireaderr check
 608
 609       END subroutine rd_grib2
 610
 611 !*****************************************************************************!
 612 ! Subroutine edition_num                                                      !
 613 !                                                                             !
 614 ! Purpose:                                                                    !
 615 !    Read one record from the input GRIB2 file.  Based on the information in  !
 616 !    the GRIB2 header and the user-defined Vtable, decide whether the field in!
 617 !    the GRIB2 record is one to process or to skip.  If the field is one we   !
 618 !    want to keep, extract the data from the GRIB2 record, and pass the data  !
 619 !    back to the calling routine.                                             !
 620 !                                                                             !
 621 ! Argument list:                                                              !
 622 !    Input:                                                                   !
 623 !       JUNIT   : "Unit Number" to open and read from.  Not really a Fortran  !
 624 !                 unit number, since we do not do Fortran I/O for the GRIB2   !
 625 !                 files.  Nor is it a UNIX File Descriptor returned from a C  !
 626 !                 OPEN statement.  It is really just an array index to the    !
 627 !                 array (IUARR) where the UNIX File Descriptor values are     !
 628 !                 stored.                                                     !
 629 !       GRIB2FILE: File name to open, if it is not already open.              !
 630 !                                                                             !
 631 !    Output:                                                                  !
 632 !       GRIB_EDITION: Set to 1 for GRIB and set to 2 for GRIB2                !
 633 !       IERR     : Error flag: 0 - no error on read from GRIB2 file.          !
 634 !                              1 - Hit the end of the GRIB2 file.             !
 635 !                              2 - The file GRIBFLNM we tried to open does    !
 636 !                                  not exist.                                 !
 637 ! Author: Paula McCaslin                                                      !
 638 ! NOAA/FSL                                                                    !
 639 ! Sept 2004                                                                   !
 640 !*****************************************************************************!
 641
 642       SUBROUTINE edition_num(junit, gribflnm,
 643      &  grib_edition, ireaderr)
 644
 645       use grib_mod
 646       use params
 647       use module_debug
 648
 649       parameter(msk1=32000,msk2=4000)
 650       character(len=1),allocatable,dimension(:) :: cgrib
 651       integer :: listsec0(3)
 652       integer :: listsec1(13)
 653       character(len=*)  :: gribflnm
 654       integer :: lskip, lgrib
 655       integer :: junit
 656       integer :: grib_edition
 657       integer :: i, j, ireaderr
 658       integer :: currlen
 659
 660       character(len=4) :: ctemp
 661       character(len=4),parameter :: grib='GRIB',c7777='7777'
 662
 663 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 664 C  SET ARGUMENTS
 665
 666       itot=0
 667       icount=0
 668       iseek=0
 669       lskip=0
 670       lgrib=0
 671       currlen=0
 672
 673 !/* IOS Return Codes from BACIO:  */
 674 !/*  0    All was well                                   */
 675 !/* -1    Tried to open read only _and_ write only       */
 676 !/* -2    Tried to read and write in the same call       */
 677 !/* -3    Internal failure in name processing            */
 678 !/* -4    Failure in opening file                        */
 679 !/* -5    Tried to read on a write-only file             */
 680 !/* -6    Failed in read to find the 'start' location    */
 681 !/* -7    Tried to write to a read only file             */
 682 !/* -8    Failed in write to find the 'start' location   */
 683 !/* -9    Error in close                                 */
 684 !/* -10   Read or wrote fewer data than requested        */
 685
 686 !if ireaderr =1 we have hit the end of a file.
 687 !if ireaderr =2 we have hit the end of all the files.
 688 !if ireaderr =3 beginning characters 'GRIB' not found
 689
 690       ! Open a byte-addressable file.
 691       CALL BAOPENR(junit,gribflnm,IOS)
 692       if (ios.eq.0) then
 693
 694          ! Search opend file for the next GRIB2 messege (record).
 695          call skgb(junit,iseek,msk1,lskip,lgrib)
 696
 697          ! Check for EOF, or problem
 698          call mprintf((lgrib.eq.0),ERROR,
 699      &     "Grib2 file or date problem, stopping in edition_num.",
 700      &     newline=.true.)
 701
 702          ! Check size, if needed allocate more memory.
 703          if (lgrib.gt.currlen) then
 704             if (allocated(cgrib)) deallocate(cgrib)
 705             allocate(cgrib(lgrib),stat=is)
 706             currlen=lgrib
 707          endif
 708
 709          ! Read a given number of bytes from unblocked file.
 710          call baread(junit,lskip,lgrib,lengrib,cgrib)
 711
 712          ! Check for beginning of GRIB message in the first 100 bytes
 713          istart=0
 714          do j=1,100
 715             ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
 716             if (ctemp.eq.grib ) then
 717               istart=j
 718               exit
 719             endif
 720          enddo
 721          if (istart.eq.0) then
 722             ireaderr=3
 723             print*, "The beginning 4 characters >GRIB< were not found."
 724          endif
 725
 726          ! Unpack Section 0 - Indicator Section to extract GRIB edition field
 727          iofst=8*(istart+5)
 728          call gbyte(cgrib,discipline,iofst,8)     ! Discipline
 729          iofst=iofst+8
 730          call gbyte(cgrib,grib_edition,iofst,8)   ! GRIB edition number
 731
 732          print *, 'ungrib - grib edition num',  grib_edition
 733          CALL BACLOSE(junit,IOS)
 734          ireaderr=1
 735       else if (ios .eq. -4) then
 736         call mprintf(.true.,ERROR,
 737      &    "edition_num: unable to open %s",newline=.true.,s1=gribflnm)
 738       else
 739          print *,'edition_num: open status failed because',ios,gribflnm
 740          ireaderr=2
 741       endif ! ireaderr check
 742
 743       END subroutine edition_num
 744
 745 !*****************************************************************************!
 746
 747       SUBROUTINE g2_compute_rh_spechumd_upa(ix, jx, iiplvl)
 748       ! Compute relative humidity from specific humidity in the upper air.
 749       use storage_module
 750       implicit none
 751       integer :: ix, jx
 752       integer :: iiplvl
 753       real :: lat1, lon1, dx, dy
 754       real, dimension(ix,jx) :: T, P, RH, Q
 755
 756       real, parameter :: svp1=611.2
 757       real, parameter :: svp2=17.67
 758       real, parameter :: svp3=29.65
 759       real, parameter :: svpt0=273.15
 760       real, parameter :: eps = 0.622
 761
 762       real startlat, startlon, deltalat, deltalon
 763
 764       call get_storage(iiplvl, 'P', P, ix, jx)
 765       call get_storage(iiplvl, 'TT', T, ix, jx)
 766       call get_storage(iiplvl, 'SH', Q, ix, jx)
 767
 768       rh=1.E2*(p*q/(q*(1.-eps)+eps))/(svp1*exp(svp2*(t-svpt0)/(T-svp3)))
 769
 770       call put_storage(iiplvl, 'RH', rh, ix, jx)
 771
 772       end subroutine g2_compute_rh_spechumd_upa
 773
 774 !*****************************************************************************!
 775
 776       subroutine ncep_grid_num (pnum)
 777 !
 778 !  Grib2 doesn't have a grid-number entry, so we have to figure it out
 779 !
 780       use gridinfo       ! Included to define map%
 781       integer :: pnum
 782       real, parameter :: eps = .01
 783       character (len=8) :: tmp8
 784
 785 !     write(6,*) 'begin ncep_grid_num'
 786 !     write(6,*) 'dx = ',map%dx,' pnum = ',pnum,' nx = ',map%nx
 787       tmp8 = '        '
 788       if (pnum .eq. 30) then
 789         if ( abs(map%dx - 12.19058) .lt. eps .and. map%nx .eq. 614) then
 790           write(tmp8,'("GRID 218")')
 791         else if (abs(map%dx - 40.63525) .lt. eps
 792      &     .and. map%nx .eq. 185) then
 793           write(tmp8,'("GRID 212")')
 794         else if (abs(map%dx - 40.63525) .lt. eps
 795      &     .and. map%nx .eq. 151) then
 796           write(tmp8,'("GRID 236")')
 797         else if (abs(map%dx - 81.2705) .lt. eps
 798      &     .and. map%nx .eq. 93) then
 799           write(tmp8,'("GRID 211")')
 800         else if (abs (map%dx - 32.46341) .lt. eps
 801      &     .and. map%nx .eq. 349) then
 802           write(tmp8,'("GRID 221")')
 803         else if (abs(map%dx - 20.317625) .lt. eps
 804      &     .and. map%nx .eq. 301) then
 805           write(tmp8,'("GRID 252")')
 806         endif
 807       else if (pnum .eq. 20) then
 808         if (abs(map%dx - 15.0) .lt. eps) then
 809           write(tmp8,'("GRID  88")')
 810         endif
 811       else if (pnum .eq. 0) then
 812         if (abs(map%dx - 1.) .lt. eps .and. map%nx .eq. 360) then
 813           write(tmp8,'("GRID   3")')
 814         else if (abs(map%dx - 0.5) .lt. eps .and. map%nx .eq. 720) then
 815           write(tmp8,'("GRID   4")')
 816         endif
 817       endif
 818       map%source(25:32) = tmp8
 819 !     write(6,*) 'map%source = ',map%source
 820       end subroutine ncep_grid_num
 821 !*****************************************************************************!
 822
 823       function earth_radius (icode)
 824 ! Grib2 Code Table 3.2. Returns the spherical earth's radius in km.
 825       use module_debug
 826       real :: earth_radius
 827       integer :: icode
 828       if ( icode .eq. 0 ) then
 829         earth_radius = 6367470. * .001
 830       else if ( icode .eq. 6 ) then
 831         earth_radius = 6371229. * .001
 832       else
 833         call mprintf(.true.,ERROR,
 834      &    "unknown earth radius for code %i",newline=.true.,i1=icode)
 835       endif
 836       end function earth_radius