!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RDNMCGRB                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE rdnmcgrb(nx_ext,ny_ext,nz_ext,                               & 8,5
           gribfile,grbflen, gribtime,                                  &
           gridesc, iproj_grb, gthin,                                   &
           ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,               &
           pi_grb,pj_grb,ipole, di_grb,dj_grb,                          &
           latsw,lonsw, latne,lonne,                                    &
           latrot,lonrot,angrot,                                        &
           latstr,lonstr,facstr,                                        &
           lattru1,lattru2,lontrue,                                     &
           scanmode, iscan,jscan,kscan,                                 &
           ires,iearth,icomp,                                           &
           jpenta,kpenta,mpenta,ispect,icoeff,                          &
           xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                         &
           rcdata,var_grb2d,var_grb3d,var_lev3d, iret)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine is to read and select given variables with GRIB
!  grid ID and level type from NMC GRIB files.
!
!  The decoder of GRIB is from NMC, ftp:nic.fb4.noaa.gov/pub
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/02/1996 Initialized.
!
!  MODIFICATIONS:
!
!  09/24/1998 (Dennis A. Moon, Chief Scientist, SSESCO)
!  Added sorting in vertical levels to assure that the level changes
!  monotonically decreasing.
!
!  1999/08/03 (Gene Bassett)
!  Excluded soil levels (grid types 111 & 112) from vertical level sorting.
!
!  2000/01/05 (Eric Kemp)
!  Fixed problem with century change over (ipds(21) runs from 1 to 100).
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    gribfile      GRIB file name
!    grbflen       Length of GRIB file name
!    gribtime      Time of GRIB data
!
!  OUTPUT:
!
!    iproj_grb     Map projection number of GRIB data
!    trlon_grb     True longitude of GRIB data (degrees E)
!    latnot_grb(2) True latitude(s) of GRIB data (degrees N)
!    swlat_grb     Latitude  of first grid point at southwest conner
!    swlon_grb     Longitude of first grid point at southwest conner
!
!    dx_grb        x-direction grid length
!    dy_grb        y-direction grid length
!
!    iret          Return flag
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INCLUDE 'gribcst.inc'      ! General GRIB definitions

  INTEGER :: nx_ext,ny_ext,nz_ext

  INTEGER :: grbflen
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  CHARACTER (LEN=1) :: opnmod
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: gdsflg       ! GDS indicator

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth

  INTEGER :: iret         ! Return flag

!
!-----------------------------------------------------------------------
!
!  Temporary arrays to read GRIB file
!
!-----------------------------------------------------------------------
!
  INTEGER :: nrecs        ! number of records

  INTEGER :: fsize        ! Size of GRIB file
  INTEGER :: grbunit      ! I/O unit of GRIB file

  LOGICAL :: ibms(nx_ext*ny_ext)    ! BMS logical array for data bit map

  REAL :: rcdata(nx_ext*ny_ext)      ! temporary data array
  REAL :: var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt) ! GRIB variables
  REAL :: var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt) ! GRIB 3-D variables
  INTEGER :: var_lev3d(nz_ext,n3dvs,n3dlvt)
                               ! Levels (hybrid) for each 3-D variable
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,l,m,n, ij, ki, kj, ilev
  INTEGER :: iyr,imo,iday,ihr,fhr, myr

  INTEGER :: ileft, iright, iincr
  INTEGER :: jleft, jright, jincr

  INTEGER :: chklev, mlvtp
  INTEGER :: iendn,itypec,wdlen, nwrd

  LOGICAL :: grdflg           ! flag for grid information
  LOGICAL :: fexist
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Find the type of endian, type of character set and length of
!  machine word
!
!  iendn     -  Integer for big-endian or little-endian
!               = 0   big-endian
!               = 1   little-endian
!               = 2   cannot compute
!  itypec    -  Integer for type of character set
!               = 0   ascii  character set
!               = 1   ebcdic character set
!               = 2   not ascii or ebcdic
!  wdlen     -  Integer for words size of computer in bytes
!               = 4   for 32 bit (4 bytes) computers
!               = 8   for 64 bit (8 bytes) computers
!
!-----------------------------------------------------------------------
!
  CALL findvartype(iendn,itypec,wdlen)
  IF ( iendn == 2 ) THEN
    WRITE(6,'(a)') 'Unknown endian type. Stopped in rdnmcgrb'
    STOP
  END IF

  IF ( itypec == 2 ) THEN
    WRITE(6,'(a)') 'Unknown character set. Stopped in rdnmcgrb'
    STOP
  END IF
!
!-----------------------------------------------------------------------
!
!   initialize the var_lev3d array so we can sort the records by
!   level values as we read them. We can't assume monotonic ordering
!   of the levels
!
!-----------------------------------------------------------------------
!
  DO k=1, nz_ext
    DO n=1, n3dvs
      DO m=1, n3dlvt
        var_lev3d(k,n,m)= -999999
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Open the GRIB file
!
!-----------------------------------------------------------------------
!
  READ (gribtime,'(i4,3i2,1x,i2)') iyr,imo,iday,ihr,fhr
  myr = MOD( iyr, 100 )

  IF (myr == 0) THEN
    myr = 100
  END IF
  INQUIRE(FILE=gribfile(1:grbflen), EXIST = fexist )
  IF ( .NOT.fexist ) THEN
    WRITE(6,'(a)') 'GRIB file '//gribfile(1:grbflen)//                  &
        ' not found. Program stopped in RDNMCGRB.'
    STOP
  END IF

  PRINT *, ' Opening GRIB file, gribfile= ',gribfile
  opnmod = 'r'
  CALL gopen  ( grbunit, gribfile, grbflen, opnmod, iret )
  IF  ( iret /= 0 )  THEN
    WRITE(6,'(a,a,/,a,i4)')                                             &
        'Error in opening file ',gribfile,                              &
        'C function gopen return status: ',iret
    RETURN
  END IF
!
!-----------------------------------------------------------------------
!
!  Get the size of the GRIB file
!
!-----------------------------------------------------------------------
!
  fsize = 0
  CALL gseek( grbunit, fsize, 2, iret )
  WRITE (6,'(a,i16,a)') 'The size of GRIB file is ',fsize,' bytes'
!
!-----------------------------------------------------------------------
!
!  Scan the GRIB file to count records
!
!-----------------------------------------------------------------------
!
  nrecs = 0
  CALL grbscan( grbunit, fsize, nprods, rcstart, rcbytes, nrecs )
  IF ( nrecs > nprods ) THEN
    WRITE (6,'(a/a,i6,a,i6/a/a)')                                       &
        'The actual number of products exceeded the expected number.',  &
        'Number of assumed products: ', nprods,                         &
        'Number of actual products: ', nrecs,                           &
        'Program stopped in RDNMCGRB. Please change NPRODS in ',         &
        'file gribcst.inc and re-run the program.'
    CALL gclose( grbunit, iret )
    STOP
  ELSE IF ( nrecs <= 0 ) THEN
    WRITE (6,'(a,a,a)')                                                 &
        'No GRIB message found in file ',gribfile,                      &
        'Program stopped in RDNMCGRB.'
    CALL gclose( grbunit, iret )
    STOP
  END IF
!
!-----------------------------------------------------------------------
!
!  Read the GRIB file and fill in the variable arrays
!
!-----------------------------------------------------------------------
!

  IF( n2dvars > n2dvs .OR. n2dlvtps > n2dlvt ) THEN
    write(6,'(a,i3,a,i3,a)')                                             &
    'Array var_nr2d(',n2dvs,',',n2dlvt,') declared in gribcst.inc too small.'
    write(6,'(a,i3,a,i3,a)')                                             &
                   'It needs to be at least ',n2dvars,'x',n2dlvtps,'.'
    write(6,'(a)') 'Program stopped in RDNMCGRB.'
    STOP
  ENDIF

  DO n=1,n2dvars
    DO m=1,n2dlvtps
      var_nr2d(n,m) = 0
    END DO
  END DO

  IF( n3dvars > n3dvs .OR. n3dlvtps > n3dlvt ) THEN
    write(6,'(a,i3,a,i3,a)')                                             &
    'Array var_nr3d(',n3dvs,',',n3dlvt,') declared in gribcst.inc too small.'
    write(6,'(a,i3,a,i3,a)')                                             &
    'It needs to be at least ',n3dvars,'x',n3dlvtps,'.'
    write(6,'(a)') 'Program stopped in RDNMCGRB.'
    STOP
  ENDIF

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_nr3d(n,m) = 0
    END DO
  END DO

  grdflg = .false.

  DO l=1,nrecs
    DO i=1,ipdsz
      ipds(i) = 0
    END DO
    DO i=1,igdsz
      igds(i) = 0
    END DO

    CALL gseek( grbunit, rcstart(l), 0, iret )
    CALL gread( grbunit, mgrib, rcbytes(l), iret )

    IF ( iendn == 1 ) THEN     ! little endian machine
      nwrd = (rcbytes(l)+3)/4
      CALL swap4byte( mgrib, nwrd )
    END IF

    CALL w3fi63( mgrib, ipds, igds, ibms, rcdata, mptrs, iret )

!    print *,'PDS#',L,ipds(3),ipds(5),ipds(6),ipds(7),rcdata(1)
    IF ( ipds(3) /= gridtyp ) THEN
!      write (6,'(a,i6,a/a,i6,a,i6)') 'Grid id of record, ',L,
!    :    'was inconsistant with expected,',
!    :    'ipds(3)=',ipds(3), ' gridtyp=',gridtyp
      CYCLE
    END IF

    DO m=1,n3dlvtps
      IF ( ipds(6) == levtyp3d(m) ) THEN
        chklev = 1
        mlvtp  = m
        GO TO 30
      ELSE
        chklev = 0
      END IF
    END DO

    DO m=1,n2dlvtps
      IF ( ipds(6) == levtyp2d(m) ) THEN
        chklev = 2
        mlvtp  = m
        EXIT
      ELSE
        chklev = 0
      END IF
    END DO

    30      CONTINUE

    IF ( chklev == 0 ) THEN
!      write (6,'(a,i4,a/a,i4)')
!    :    'Type of level of record, ',L,
!    :    ' was inconsistant with expected,',
!    :    ' ipds(6) =',ipds(6)
      CYCLE
    END IF

    IF ( ipds(8) /= myr .OR.ipds(9) /= imo.OR.                          &
           ipds(10) /= iday.OR.ipds(11) /= ihr.OR.                      &
           ipds(14) /= fhr ) THEN
      WRITE (6,'(a,i6,a/a,4i2.2,i2.2/a)')                               &
          'Time of record, ',l,                                         &
          'was inconsistant with expected,',                            &
          ' pds time = ',(ipds(j),j=8,11),ipds(14),                     &
          ' gribtime = ',gribtime
      CYCLE
    END IF

    gdsflg = IAND(ipds(4),128)/128
    IF ( .NOT. grdflg ) THEN
      CALL grbgrid(gridtyp, gdsflg, igds,                               &
                   gridesc, iproj_grb, gthin,                           &
                   ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,       &
                   pi_grb,pj_grb,ipole, di_grb,dj_grb,                  &
                   latsw,lonsw, latne,lonne,                            &
                   latrot,lonrot,angrot,                                &
                   latstr,lonstr,facstr,                                &
                   lattru1,lattru2,lontrue,                             &
                   scanmode, iscan,jscan,kscan,                         &
                   ires,iearth,icomp,                                   &
                   jpenta,kpenta,mpenta,ispect,icoeff,                  &
                   xp_grb,yp_grb, xo_grb,yo_grb,zo_grb)

      WRITE (6,'(a)') gridesc

      grdflg = .true.
    END IF

    IF ( igds(1) /= mproj_grb ) THEN
!      write (6,'(a,i3)')
!    :    'Map projection was not consistant, igds(1) = ', igds(1)
      CYCLE
    END IF

    IF ( igds(2) /= nx_ext .OR. igds(3) /= ny_ext) THEN
      WRITE (6,'(a/a,i4,2x,a,i4/a,i4,2x,a,i4/a)')                       &
          'Horizontal dimension was not consistant,',                   &
          ' igds(2) = ', igds(2), ' igds(3) = ', igds(3),               &
          ' nx_ext  = ', nx_ext  , ' ny_ext  = ', ny_ext,               &
          'Program stopped in RDNMCGRB.'
      STOP
    END IF

    IF ( iscan == 0 ) THEN
      ileft  = 1
      iright = nx_ext
      iincr  = 1
    ELSE
      ileft  = nx_ext
      iright = 1
      iincr  = -1
    END IF

    IF ( jscan /= 0 ) THEN
      jleft  = 1
      jright = ny_ext
      jincr  = 1
    ELSE
      jleft  = ny_ext
      jright = 1
      jincr  = -1
    END IF

    IF ( kscan == 0 ) THEN
      ki = 1
      kj = nx_ext
    ELSE IF ( kscan == 1 ) THEN
      ki = ny_ext
      kj = 1
    ELSE
!      write (6,'(a,i3,a,i6)')
!    :    'Unknown scanning mode, kScan = ', kScan
      CYCLE
    END IF

    IF ( chklev == 1 ) THEN           ! 3-D variable
      DO n=1,n3dvars
        DO m=1,n3dlvtps

          IF ( ipds(5) == var_id3d(n,m) .AND.  ipds(6) == levtyp3d(m) ) THEN
            var_nr3d(n,m) = var_nr3d(n,m) + 1

!          print '(5(a,i6))',
!    :        'ipds(7) = ', ipds(7),' var = ',var_id3d(n,m),
!    :        ' var_nr3d = ',var_nr3d(n,m),' n = ', n, ' m = ', m
!

            IF (levtyp3d(m) == 111 .OR. levtyp3d(m) == 112) THEN
!
!-----------------------------------------------------------------------
!
!  Insert the level and data values into the var_lev3d and var_grb3d
!  in the order that they occur in the data file (e.g. below surface
!  coordinates).
!
!-----------------------------------------------------------------------
!
              ilev = var_nr3d(n,m)

            ELSE
!
!-----------------------------------------------------------------------
!
!  Insert the level and data values into the var_lev3d and var_grb3d
!  arrays in order of decreasing level value,
!
!  For example if the level type is pressure levels the levels
!  will be sorted in order of decreasing pressure
!
!  ilev is the level at which to insert the current data
!
!-----------------------------------------------------------------------
!
              ilev= 0
              DO k=1,nz_ext
                IF ( ipds(7) > var_lev3d(k,n,m) ) THEN
                  ilev = k
                  EXIT
                END IF
              END DO

!              120             CONTINUE

              IF ( ilev == 0 ) THEN
                PRINT*,'couldnt locate level ',ipds(7)
                PRINT*,'for var ID#',var_id3d(n,m)
                PRINT *, ' n = ', n, ' m = ', m
                STOP
              END IF

!   slide all the following levels down one to make room

              DO k=nz_ext,ilev+1, -1
                var_lev3d(k,n,m)= var_lev3d(k-1,n,m)
                DO j=jleft,jright,jincr
                  DO i=ileft,iright,iincr
                    var_grb3d(i,j,k,n,m)= var_grb3d(i,j,k-1,n,m)
                  END DO
                END DO
              END DO

            END IF

!
!   now insert the data at the appropriate level
!

            var_lev3d(ilev,n,m) = ipds(7)
            DO j=jleft,jright,jincr
              DO i=ileft,iright,iincr
                ij = ki*iincr*(i-ileft) + kj*jincr*(j-jleft) + 1
                var_grb3d(i,j,ilev,n,m) = rcdata(ij)
              END DO
            END DO
          END IF
        END DO
      END DO
    ELSE IF ( chklev == 2 ) THEN       ! 2-D variable
      IF(ipds(16) == 0 .OR. ipds(16) == 10) THEN
                                       ! Don't use average, accumulation,
! or time difference fields.
        DO n=1,n2dvars
          DO m=1,n2dlvtps
            IF ( ipds(5) == var_id2d(n,m) ) THEN
              var_nr2d(n,m) = var_nr2d(n,m) + 1
              DO i=ileft,iright,iincr
                DO j=jleft,jright,jincr
                  ij = ki*iincr*(i-ileft) + kj*jincr*(j-jleft) + 1
                  var_grb2d(i,j,n,m) = rcdata(ij)
                END DO
              END DO
            END IF
          END DO
        END DO
      END IF
    END IF

  END DO

!
!-----------------------------------------------------------------------
!
!  Set good status
!
!-----------------------------------------------------------------------
!
  iret = 0

  CALL gclose( grbunit, iret )

  RETURN
END SUBROUTINE rdnmcgrb
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GRBSCAN                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE grbscan(grbunit,fsize,nprods,                                & 1
           rcstart,rcbytes,nrecs)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine is to scan the NMC GRIB files.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  01/03/1996 Initialized.
!
!  MODIFICATIONS:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    grbunit  Unit of GRIB file opened to scan for GRIB records
!    nprods   Maxmum record number
!    fsize    Size of the GRIB file
!
!  OUTPUT:
!
!    rcstart  Starting position (bytes) for each record
!    rcbytes  Record length in bytes
!    nrecs    Record number actually scanned
!    iret     Return flag
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: grbunit
  INTEGER :: fsize
  INTEGER :: nprods

  INTEGER :: rcbytes(nprods)  ! record length in bytes
  INTEGER :: rcstart(nprods)  ! record starting byte in a GRIB file

  INTEGER :: nrecs
  INTEGER :: iret
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nbytes

  CHARACTER (LEN=1) :: buffer(8)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nbytes = 0
  nrecs  = 0

  100   IF ( nbytes < fsize ) THEN
    CALL gseek( grbunit,nbytes,0,iret )
    IF ( iret /= 0 ) GO TO 200

    CALL gread( grbunit,buffer,8,iret )
    IF ( iret /= 8 ) GO TO 200

    IF ( ICHAR(buffer(1)) /= 71 .OR.     & ! test ASCII 'G'
     ICHAR(buffer(2)) /= 82 .OR.     & ! test ASCII 'R'
     ICHAR(buffer(3)) /= 73 .OR.     & ! test ASCII 'I'
     ICHAR(buffer(4)) /= 66 .OR.     & ! test ASCII 'B'
     ICHAR(buffer(8)) /= 1 ) THEN   ! test Edition #, 1
    nbytes = nbytes + 1
    GO TO 100
  END IF

  nrecs = nrecs + 1
  IF ( nrecs <= nprods ) THEN
    rcstart(nrecs) = nbytes
    rcbytes(nrecs) = ICHAR(buffer(5))*65536                             &
                   + ICHAR(buffer(6))*256                               &
                   + ICHAR(buffer(7))
    nbytes = nbytes + rcbytes(nrecs)
  END IF

  GO TO 100

  END IF

  200   RETURN

END SUBROUTINE grbscan