/WPS/ungrib/src/gribcode.F

!                                                                             !
!*****************************************************************************!
!                                                                             !
! This is a package of subroutines to read GRIB-formatted data.  It is still  !
! under continuous development. It will not be able to read every GRIB dataset!
! you give it, but it will read a good many.                                  !
!                                                                             !
!   Kevin W. Manning                                                          !
!   NCAR/MMM                                                                  !
!   Summer 1998, and continuing                                               !
!   SDG                                                                       !
!                                                                             !
!*****************************************************************************!
!                                                                             !
! The main user interfaces are:                                               !
!                                                                             !
!  SUBROUTINE GRIBGET(NUNIT, IERR)                                            !
!      Read a single GRIB record from UNIX file-descriptor NUNIT into array   !
!      GREC. No unpacking of any header or data values is performed.          !
!                                                                             !
!  SUBROUTINE GRIBREAD(NUNIT, DATA, NDATA, IERR)                              !
!      Read a single GRIB record from UNIX file-descriptor NUNIT, and unpack  !
!      all header and data values into the appropriate arrays.                !
!                                                                             !
!  SUBROUTINE GRIBHEADER                                                      !
!      Unpack the header of a GRIB record                                     !
!                                                                             !
!  SUBROUTINE GRIBDATA(DATARRAY, NDAT)                                        !
!      Unpack the data in a GRIB record into array DATARRAY                   !
!                                                                             !
!  SUBROUTINE GRIBPRINT(ISEC)                                                 !
!      Print the header information from GRIB section ISEC.                   !
!                                                                             !
!  SUBROUTINE GET_SEC1(KSEC1)                                                 !
!      Return the header information from Section 1.                          !
!                                                                             !
!  SUBROUTINE GET_SEC2(KSEC2)                                                 !
!      Return the header information from Section 2.                          !
!                                                                             !
!  SUBROUTINE GET_GRIDINFO(IGINFO, GINFO)                                     !
!      Return the grid information of the previously-unpacked GRIB header.    !
!                                                                             !
!                                                                             !
!*****************************************************************************!
!                                                                             !
!                                                                             !
! The following arrays have meanings as follows:                              !
!                                                                             !
!                                                                             !
! SEC0:    GRIB Header Section 0 information                                  !
!                                                                             !
!       1  :  Length of a complete GRIB record                                !
!       2  :  GRIB Edition number                                             !
!                                                                             !
!                                                                             !
! SEC1:    GRIB Header Section 1 information                                  !
!                                                                             !
!       1  :  Length of GRIB section 1 (bytes)                                !
!       2  :  Parameter Table Version number ????                             !
!       3  :  Center Identifier ????                                          !
!       4  :  Process Identifier ????                                         !
!       5  :  Grid ID number for pre-specified grids.                         !
!       6  :  Binary bitmap flag:                                             !
!       7  :  Parameter ID Number (ON388 Table 2)                             !
!       8  :  Level type (ON388 Table 3)                                      !
!       9  :  Level value, or top value of a layer                            !
!      10  :  Bottom value of a layer ( 0 if NA ??)                           !
!      11  :  Year (00-99)                                                    !
!      12  :  Month (01-12)                                                   !
!      13  :  Day of the month (01-31)                                        !
!      14  :  Hour (00-23)                                                    !
!      15  :  Minute (00-59)                                                  !
!      16  :  Forecast time unit: (ON388 Table 4)                             !
!      17  :  Time period 1:                                                  !
!      18  :  Time period 2:                                                  !
!      19  :  Time range indicator (ON833 Table 5)                            !
!      20  :  Number of ?? in an average ??                                   !
!      21  :  Number of ?? missing from average ??                            !
!      22  :  Century (Years 1999 and 2000 are century 20, 2001 is century 21)!
!      23  :  Sub-center identifier ??                                        !
!      24  :  Decimal scale factor for ???                                    !
!                                                                             !
!                                                                             !
!                                                                             !
!                                                                             !
!                                                                             !
! SEC2:    GRIB Header Section 2 information                                  !
!                                                                             !
!       1  :  Length of GRIB Section 2                                        !
!       2  :  Number of vertical-coordinate parameters ???                    !
!       3  :  Starting-point of the list of vertical-coordinate parameters ?? !
!       4  :  Data-representation type (i.e., grid type)  Table ???           !
!          :      0 = ??                                                      !
!          :      3 = Lambert-conformal grid.                                 !
!          :      5 = Polar-stereographic grid.                               !
!                                                                             !
!     if (SEC2(4) == 0) then            LATITUDE/LONGITUDE GRID               !
!                                                                             !
!         INFOGRID/GRIDINFO:                                                  !
!                                                                             !
!             1  :  I Dimension of the grid                                   !
!             2  :  J Dimension of the grid                                   !
!             3  :  Starting Latitude of the grid.                            !
!             4  :  Starting Longitude of the grid.                           !
!             5  :  Resolution and component flags.                           !
!             6  :  Ending latitude of the grid.                              !
!             7  :  Ending longitude of the grid.                             !
!             8  :  Longitudinal increment.                                   !
!             9  :  Latitudinal incriment.                                    !
!            10  :  Scanning mode (bit 3 from Table 8)                        !
!            21  :  Iscan sign (+1/-1) (bit 1 from Table 8)                   !
!            22  :  Jscan sign (+1/-1) (bit 2 from Table 8)                   !
!                                                                             !
!                                                                             !
!     elseif (SEC2(4) == 3) then        LAMBERT CONFORMAL GRID                !
!                                                                             !
!         INFOGRID/GRIDINFO:                                                  !
!                                                                             !
!             1  :  I Dimension of the grid                                   !
!             2  :  J Dimension of the grid                                   !
!             3  :  Starting Latitude of the grid.                            !
!             4  :  Starting Longitude of the grid.                           !
!             5  :  Resolution and component flags.                           !
!             6  :  Center longitude of the projection.                       !
!             7  :  Grid-spacing in the I direction                           !
!             8  :  Grid-spacing in the J direction                           !
!             9  :  Projection center                                         !
!            10  :  Scanning mode (bit 3 from Table 8)                        !
!            11  :  First TRUELAT value.                                      !
!            12  :  Second TRUELAT value.                                     !
!            13  :  Latitude of the southern pole ??                          !
!            14  :  Longitude of the southern pole ??                         !
!            21  :  Iscan sign (+1/-1) (bit 1 from Table 8)                   !
!            22  :  Jscan sign (+1/-1) (bit 2 from Table 8)                   !
!                                                                             !
!     if (SEC2(4) == 4) then            GAUSSIAN GRID                         !
!                                                                             !
!         INFOGRID/GRIDINFO:                                                  !
!                                                                             !
!             1  :  I Dimension of the grid                                   !
!             2  :  J Dimension of the grid                                   !
!             3  :  Starting Latitude of the grid.                            !
!             4  :  Starting Longitude of the grid.                           !
!             5  :  Resolution and component flags.                           !
!             6  :  Ending latitude of the grid.                              !
!             7  :  Ending longitude of the grid.                             !
!             8  :  Longitudinal increment.                                   !
!             9  :  Number of latitude circles between pole and equator       !
!            10  :  Scanning mode (bit 3 from Table 8)                        !
!            17  :  Original (stored) ending latitude                         !
!            18  :  Original (stored) starting latitude                       !
!            19  :  Approximate delta-latitude                                !
!            21  :  Iscan sign (+1/-1) (bit 1 from Table 8)                   !
!            22  :  Jscan sign (+1/-1) (bit 2 from Table 8)                   !
!                                                                             !
!                                                                             !
!     elseif (SEC2(4) == 5) then        POLAR STEREOGRAPHIC GRID              !
!                                                                             !
!         INFOGRID/GRIDINFO:                                                  !
!                                                                             !
!             1  :  I Dimension of the grid                                   !
!             2  :  J Dimension of the grid                                   !
!             3  :  Starting Latitude of the grid.                            !
!             4  :  Starting Longitude of the grid.                           !
!             5  :  Resolution and component flags.                           !
!             6  :  Center longitude of the projection.                       !
!             7  :  Grid-spacing in the I direction                           !
!             8  :  Grid-spacing in the J direction                           !
!             9  :  Projection center                                         !
!            10  :  Scanning mode (bit 3 from Table 8)                        !
!            21  :  Iscan sign (+1/-1) (bit 1 from Table 8)                   !
!            22  :  Jscan sign (+1/-1) (bit 2 from Table 8)                   !
!                                                                             !
!     elseif (SEC2(4) == 50) then       SPHERICAL HARMONIC COEFFICIENTS       !
!                                                                             !
!         INFOGRID/GRIDINFO:                                                  !
!                                                                             !
!             1  :  J-pentagonal resolution parameter                         !
!             2  :  K-pentagonal resolution parameter                         !
!             3  :  M-pentagonal resolution parameter                         !
!             4  :  Spectral representation type (ON388 Table 9)              !
!             5  :  Coefficient storage mode (ON388 Table 10)                 !
!                                                                             !
!     elseif (SEC2(4) == ?) then        ??                                    !
!                                                                             !
!                                                                             !
! SEC3:    GRIB Header Section 3 information:                                 !
! SEC4:    GRIB Header Section 4 information:                                 !
!
!
!
module module_grib
!
! Machine wordsize must be known for the various unpacking routines to work.
! Machine wordsize is set through CPP Directives.
! Use options -DBIT32 (for 32-bit word-size) or -DBIT64 (for 64-bit wordsize) 
! for the CPP pass of the compiler.
!
#if defined (BIT32)
  integer, parameter :: MWSIZE = 32 ! Machine word size in bits
#elif defined (BIT64)
  integer, parameter :: MWSIZE = 64 ! Machine word size in bits
#endif


! Array GREC holds a single packed GRIB record (header and all).
! Array BITMAP holds the bitmap (if a bitmap was used).
!
! For some reason, the cray does not like grec to be allocatable.
!
#if defined (CRAY) 
  integer, dimension(100000) :: grec
  integer, dimension(100000) :: bitmap
#else
  integer, allocatable, save, dimension(:) :: grec
  integer, allocatable, save, dimension(:) :: bitmap
#endif

! SEC0 holds the Section 0 header information
! SEC1 holds the Section 1 header information
! SEC2 holds the Section 2 header information
! SEC3 holds the Section 3 header information
! SEC4 holds the Section 4 header information
! XEC4 holds floating-point Section 4 header information

  integer, dimension(2) :: sec0
  integer, dimension(100) :: sec1
  integer, dimension(10) :: sec2
  integer, dimension(10) :: sec3
  integer, dimension(10) :: sec4
  real, dimension(1) :: xec4

  integer :: sevencount = 0

! INFOGRID holds integer values defining the grid.
! GRIDINFO holds floating-point values definint the grid

  integer, dimension(40) :: infogrid
  real, dimension(40) :: gridinfo

  integer :: ied
  real, parameter :: pi = 3.1415926534
  real, parameter :: degran = pi/180.
  real, parameter :: raddeg = 1./degran

  real :: glat1, glon1, gclon, gtrue1, gtrue2, grrth, gx1, gy1, gkappa

contains
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
  integer function gribsize(trec, ilen, ierr)
!-----------------------------------------------------------------------------!
! Return the size of a single GRIB record.                                    !
!                                                                             !
! Input:                                                                      !
!    TREC: At least a portion of the complete GRIB record.                    !
!    ILEN: The size of array TREC.                                            !
!                                                                             !
! Output                                                                      !
!    GRIBSIZE: The size of the full GRIB record                               !
!    IERR    : 0= no errors, 1 = read error
!                                                                             !
! Side Effects:                                                               !
!    * Module variable IED is set to the GRIB Edition number.                 !
!    * STOP, if not GRIB Edition 0 or 1                                      !
!                                                                             !
! Externals                                                                   !
!    GBYTE                                                                    !
!                                                                             !
!-----------------------------------------------------------------------------!
    implicit none
    integer :: ilen
    integer, dimension(ilen) :: trec
    integer :: isz0 = 32
    integer :: isz1 = 0
    integer :: isz2 = 0
    integer :: isz3 = 0
    integer :: isz4 = 0
    integer :: isz5 = 32
    integer :: iflag
    integer :: ierr
    character :: pname*132

    ierr = 0
! Unpack the GRIB Edition number, located in the eighth byte (bits 57-64)     !
! of array TREC.                                                              !

    call gbyte_g1(trec, ied, 56, 8)

! GRIB Edition 1 has the size of the whole GRIB record right up front.        !

    if (ied.eq.1) then
       ! Grib1
       ! Find the size of the whole GRIB record
       call gbyte_g1(trec, gribsize, 32, 24)

! GRIB Edition 0 does not include the total size, so we have to sum up        !
! the sizes of the individual sections                                        !

    elseif (ied.eq.0) then
       ! Grib old edition
       ! Find the size of section 1.
       call gbyte_g1(trec, isz1, isz0, 24)
       isz1 = isz1 * 8
       call gbyte_g1(trec, iflag, isz0+56, 8)
       if ((iflag.eq.128).or.(iflag.eq.192)) then ! section 2 is there
          ! Find the size of section 2.
          call gbyte_g1(trec, isz2, isz0+isz1, 24)
          isz2 = isz2 * 8
       endif
       if ((iflag.eq.64).or.(iflag.eq.192)) then ! Section 3 is there
          ! Find the size of section 3.
          call gbyte_g1(trec, isz3, isz0+isz1+isz2, 24)
          isz3 = isz3 * 8
       endif
       ! Find the size of section 4.
       call gbyte_g1(trec, isz4, isz0+isz1+isz2+isz3, 24)
       isz4 = isz4 * 8

       ! Total the sizes of sections 0 through 5.
       gribsize = (isz0+isz1+isz2+isz3+isz4+isz5) / 8

    elseif (ied.eq.2) then
       ! Grib2
       CALL getarg ( 0 , pname )
       write(*,'("*** stopping in gribcode ***\n")')
       write(*,'("\tI was expecting a Grib1 file, but this is a Grib2 file.")')
       if ( index(pname,'ungrib.exe') .ne. 0 ) then
         write(*,'("\tIt is possible this is because your GRIBFILE.XXX files")')
         write(*,'("\tare not all of the same type.")')
         write(*,'("\tWPS can handle both file types, but a separate ungrib")')
         write(*,'("\tjob must be run for each Grib type.\n")')
       else
         write(*,'("\tUse g2print on Grib2 files\n")')
       endif
       stop 'gribsize in gribcode'
    else
       write(*,'("Error trying to read grib edition number in gribsize.")')
       write(*,'("Possible corrupt grib file.")')
       write(6,*) 'Incorrect edition number  = ',ied
       write(6,*) 'Skipping the rest of the file and continuing.'
       ierr = 1
    endif
  end function gribsize
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
  subroutine findgrib(nunit, isize, ierr)

!-----------------------------------------------------------------------------!
!                                                                             !
! Find the string "GRIB", which starts off a GRIB record.                     !
!                                                                             !
! Input:                                                                      !
!    NUNIT:  The C unit to read from.  This should already be opened.         !
!                                                                             !
! Output:                                                                     !
!    ISIZE:  The size in bytes of one complete GRIB Record                    !
!    IERR:   Error flag,                                                      !
!              0 : No error or end-of-file on reading                         !
!              1 : Hit the end of file                                        !
!              2 : Error on reading                                           !
!                                                                             !
! Side effects:                                                               !
!   * The pointer to C unit NUNIT is set to the beginning of the next         !
!     GRIB record.                                                            !
!   * The first time FINDGRIB is called, the integer GTEST is set to          !
!     a value equivalent to the string 'GRIB'                                 !
!                                                                             !
! Modules:                                                                    !
!     MODULE_GRIB                                                             !
!                                                                             !
! Externals:                                                                  !
!     BN_READ                                                                 !
!     BN_SEEK                                                                 !
!     GRIBSIZE                                                                !
!                                                                             !
!-----------------------------------------------------------------------------!
    implicit none
    integer, intent(in) :: nunit
    integer, intent(out) :: isize
    integer, intent(out) :: ierr

    integer, parameter :: LENTMP=100
    integer, dimension(lentmp) :: trec

    integer :: isz, itest, icnt

    integer, save :: gtest = 0

! Set the integer variable GTEST to hold the integer representation of the
! character string 'GRIB'.   This integer variable is later compared to
! integers we read from the GRIB file, to find the beginning of a GRIB record.

    if (gtest.eq.0) then
       if (mwsize.eq.32) then
          gtest = transfer('GRIB', gtest)
       elseif(mwsize.eq.64) then
          call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'GRIB', gtest, 0, mwsize)
       endif
    endif
    ierr = 0
    icnt = 0

    LOOP : DO
! Read LENTMP bytes into holding array TREC.
       call bn_read(nunit, trec, lentmp, isz, ierr, 0)
       if (ierr.eq.1) then
          return
       elseif (ierr.eq.2) then
          write(*,'("Error reading GRIB: IERR = ", I2)') ierr
          return
       endif
! Reposition the file pointer back to where we started.
       call bn_seek(nunit, -isz, 0, 0)

! Compare the first four bytes of TREC with the string 'GRIB' stored in 
! integer variable GTEST.
       if (mwsize.eq.32) then
          if (trec(1) == gtest) exit LOOP
       elseif (mwsize.eq.64) then
          call gbyte_g1(trec, itest, 0, 32)
          if (itest == gtest) exit LOOP
       endif

! Advance the file pointer one byte.
       call bn_seek(nunit, 1, 0, 0)
       icnt = icnt + 1
       if ( icnt .gt. 100000) then       ! stop if we cannot find the GRIB string
         write(*,'("*** stopping in findgrib in gribcode ***\n")')
         write(*,'("\tI could not find the GRIB string in the input file")')
         write(*,'("\tafter testing the first 100,000 bytes.")')
         write(*,'("\tEither the file is corrupt or it is not a GRIB file.\n")')
         stop 'findgrib'
       endif

    ENDDO LOOP

!#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
#ifdef BYTESWAP
      call swap4(trec, isz)
#endif
    isize = gribsize(trec, isz, ierr)

  end subroutine findgrib
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
  subroutine SGUP_NOBITMAP(datarray, ndat)
! Simple grid-point unpacking
    implicit none

    integer :: ndat
    real , dimension(ndat) :: datarray
    integer, dimension(ndat) :: IX
    real :: dfac, bfac
    integer :: iskip

    DFAC = 10.**(-sec1(24))
    BFAC = 2.**sec4(7)
    if (ied.eq.0) then
       iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    elseif (ied.eq.1) then
       iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    endif
! sec4(8) is the number of bits used per datum value.
! If sec4(8) = 255, assume they mean sec4(8) = 0
    if (sec4(8) == 255) then
       sec4(8) = 0
    endif
! If sec4(8) is 0, assume datarray is constant value of xec4(1)

    if (sec4(8).eq.0) then
       !!! HERE IS THE ORIGINAL NCAR CODE: 
       ! datarray = xec4(1)
       !!! HERE IS WHAT FSL CHANGED IT TO:
       datarray = DFAC*xec4(1)
       !!! because even though it is a constant value
       !!! your still need to scale by the decimal scale factor.
    else
       !!! FSL developers MOVED THE CALL TO gbytes FROM line 441 ABOVE 
       !!! BECAUSE IF sec4(8)==0 BEFORE gbytes IS CALLED, THE MASKS ARRAY
       !!! IN gbytes WILL BE INDEXED OUT OF BOUNDS. C HARROP 9/16/04
       call gbytes_g1(grec, IX, iskip, sec4(8), 0, ndat)
       datarray = DFAC * (xec4(1) + (IX*BFAC))
    endif
  end subroutine SGUP_NOBITMAP
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!

  subroutine SGUP_BITMAP(datarray, ndat)
! Simple grid-point unpacking, with a bitmap.
    implicit none

    integer :: ndat ! Number of data points in the final grid.
    real , dimension(ndat) :: datarray ! Array holding the final unpacked data.
    real :: dfac, bfac
    integer :: iskip, nbm, i, nn

    integer, allocatable, dimension(:) :: bmdat

! SEC4(1) : The number of bytes in the whole of GRIB Section 4.
! SEC4(6) : The number of unused bits at the end of GRIB Section 4.
! SEC4(8) : The number of bits per data value.

    datarray = -1.E30

! 1) There are fewer than NDAT data values, because a bitmap was used.  
!    Compute the number of data values (NBM).  There are 11 extra bytes
!    in the header section 4.  NBM equals the total number of data bits (not
!    counting the header bits), minus the number of unused buts, and then
!    divided by the number of bits per data value.

! If sec4(8) is 0, assume datarray is constant value of xec4(1)

    if (sec4(8).eq.0) then
       where(bitmap(1:ndat).eq.1) datarray = xec4(1)
       return
    endif
    nbm = ((sec4(1)-11)*8-sec4(6))/sec4(8)
    allocate(bmdat(nbm))

! Compute the parameters involved with packing
    DFAC = 10.**(-sec1(24))
    BFAC = 2.**sec4(7)

! Set ISKIP to the beginning of the data.
    if (ied.eq.0) then
       iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    elseif (ied.eq.1) then
       iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    endif

! Read the data from the GREC array
    call gbytes_g1(grec, bmdat, iskip, sec4(8), 0, nbm)
! sec4(8) is the number of bits used per datum value.
! If sec4(8) = 255, assume they mean sec4(8) = 0
    if (sec4(8) == 255) sec4(8) = 0

! Unpack the data according to packing parameters DFAC, BFAC, and XEC4(1), 
! and masked by the bitmap BITMAP.
       nn = 0
       do i = 1, ndat
          if (bitmap(i).eq.1) then
             nn = nn + 1
             datarray(i) = DFAC * (xec4(1) + (bmdat(nn)*BFAC))
          endif
       enddo

! Deallocate the scratch BMDAT array
    deallocate(bmdat)

  end subroutine SGUP_BITMAP
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!

  subroutine CSHUP(pdata, ndat)
! ECMWFs unpacking of ECMWFs Complex Spherical Harmonic packing
! Adapted from ECMWFs GRIBEX package.
    implicit none

    integer :: ndat
    real , dimension(ndat) :: pdata
    integer, dimension(ndat+500) :: IX

    integer :: iskip, isign
    integer :: N1, IPOWER, J, K, M, nval
    real :: zscale, zref

    integer :: ic, jm, iuc, il2, inum, jn
    integer :: inext, ilast, ioff, jrow, index, i, jcol
    real :: bval
    integer, allocatable, dimension(:) :: iexp, imant
    real , dimension(0:400) :: factor
    real :: power
    integer :: N

    index = -1

    if (ied.eq.0) then
       iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    elseif(ied.eq.1) then
       iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    endif

    call gbyte_g1(grec,N1,iskip,16)
    iskip = iskip + 16

    call gbyte_g1(grec,ipower,iskip,16)
    iskip = iskip + 16
    if (ipower.ge.32768) ipower = 32768-ipower

! Unpack the resolution parameters for the initial (small) truncation:
    call gbyte_g1(grec,J,iskip,8)
    iskip = iskip + 8
    call gbyte_g1(grec,K,iskip,8)
    iskip = iskip + 8
    call gbyte_g1(grec,M,iskip,8)
    iskip = iskip + 8

    zscale = 2.**sec4(7)

    iskip = N1*8

    nval = NDAT - (J+1)*(J+2)

    call gbytes_g1(grec, ix, iskip, sec4(8), 0, nval)
! sec4(8) is the number of bits used per datum value.
! If sec4(8) = 255, assume they mean sec4(8) = 0
    if (sec4(8) == 255) sec4(8) = 0

    pdata(1:nval) = (float(ix(1:nval))*zscale)+xec4(1)

    IUC = NDAT+1
    IC  = NVAL+1
    DO JM=INFOGRID(1),0,-1
       IL2=MAX(JM,J+1)
       INUM=2*(INFOGRID(1)-IL2+1)
       pdata(iuc-inum:iuc-1) = pdata(ic-inum:ic-1)
       iuc = iuc - inum
       ic = ic - inum
       IUC = IUC-MAX((IL2-JM)*2,0)
    ENDDO

    if (ied.eq.0) then
       iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 11*8
    elseif (ied.eq.1) then
       iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8 + 18*8
    endif

    allocate(iexp(802))
    allocate(imant(802))
    ilast=j+1
    do jrow=1,ilast
       inext = 2*(ilast-jrow+1)
       ! extract all the exponents
       call gbytes_g1(grec, iexp, iskip, 8, 24, inext)
       ! extract all the mantissas
       ioff = 8
       call gbytes_g1(grec, imant, iskip+8, 24, 8, inext)
       iskip = iskip + inext*32

       ! Build the real values from mantissas and exponents
       bval = 2.**(-24)
       i = 0
       do jcol = jrow, infogrid(1)+1
          index = index + 2
          if (ilast.ge.jcol) then
             i = i + 1
             if ((iexp(i).eq.128.or.iexp(i).eq.0).and.(imant(i).eq.0)) then
                pdata(i) = 0
             else
                if (iexp(i).ge.128) then
                   iexp(i) = iexp(i) - 128
                   isign = -1
                else
                   isign = 1
                endif
                pdata(index) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
                i = i + 1
                if (iexp(i).ge.128) then
                   iexp(i) = iexp(i) - 128
                   isign = -1
                else
                   isign = 1
                endif
                pdata(index+1) = isign*bval*IMANT(i)*16.**(IEXP(i)-64)
             endif
          endif
       enddo
    enddo

    !Apply power scaling:

    if (ipower.ne.0) then
       power = float(ipower) / 1000.0
       factor(0) = 1.0
       do n = 1 , infogrid(1)
          if( ipower .ne. 1000 ) then
             factor(n) = 1.0 / (n * (n+1) )**power
          else
             factor(n) = 1.0 / (n * (n + 1))
          endif
       enddo
       INDEX = -1
       DO M = 0 , J-1
          DO N = M , INFOGRID(1)
             INDEX = INDEX + 2
             IF ( N .GE. J ) THEN
                PDATA(INDEX:INDEX+1) = PDATA(INDEX:INDEX+1) * FACTOR(N)
             ENDIF
          ENDDO
       ENDDO
       DO M = J , INFOGRID(1)
          DO N = M , INFOGRID(1)
             INDEX = INDEX + 2
             PDATA(INDEX:INDEX+1)   = PDATA(INDEX:INDEX+1)   * FACTOR(N)
          ENDDO
       ENDDO
    endif

  end subroutine CSHUP
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
!
!
! Trigonometric functions which deal with degrees, rather than radians:
!
  real function sind(theta)
    real :: theta
    sind = sin(theta*degran)
  end function sind
  real function cosd(theta)
    real :: theta
    cosd = cos(theta*degran)
  end function cosd
  real function tand(theta)
    real :: theta
    tand = tan(theta*degran)
  end function tand
  real function atand(x)
    real :: x
    atand = atan(x)*raddeg
  end function atand
  real function atan2d(x,y)
    real :: x,y
    atan2d = atan2(x,y)*raddeg
  end function atan2d
  real function asind(x)
    real :: x
    asind = asin(x)*raddeg
  end function asind
  real function acosd(x)
    real :: x
    acosd = acos(x)*raddeg
  end function acosd

!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
end module module_grib
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribget(nunit, ierr)
  use module_grib
!-----------------------------------------------------------------------------!
!                                                                             !
! Read a single GRIB record, with no unpacking of any header or data fields.  !
!                                                                             !
! Input:                                                                      !
!     NUNIT:  C unit number to read from.  This should already be open.       !
!                                                                             !
! Output:                                                                     !
!     IERR: Error flag, Non-zero means there was a problem with the read.     !
!                                                                             !
! Side Effects:                                                               !
!        The array GREC is allocated, and filled with one GRIB record.        !
!        The C unit pointer is moved to the end of the GRIB record just read. !
!                                                                             !
! Modules:                                                                    !
!       MODULE_GRIB                                                           !
!                                                                             !
! Externals:                                                                  !
!       FINDGRIB                                                              !
!       BN_READ                                                               !
!                                                                             !
!-----------------------------------------------------------------------------!

  implicit none

  integer :: nunit
  integer :: ierr
  integer :: isz, isize

! Position the file pointer at the beginning of the GRIB record.
  call findgrib(nunit, isize, ierr)
  if (ierr.ne.0) return

! Allocate the GREC array to be able to hold the data

#if defined (CRAY)
#else
  allocate(grec((isize+(mwsize/8-1))/(mwsize/8)))
#endif

! Read the full GRIB record.

  call bn_read(nunit, grec, isize, isz, ierr, 1)

!#if defined (DEC) || defined (ALPHA) || defined (alpha) || defined (LINUX)
#ifdef BYTESWAP
      call swap4(grec, isz)
#endif


end subroutine gribget
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribread(nunit, data, ndata, debug_level, ierr)
!-----------------------------------------------------------------------------!
! Read one grib record, unpack the header and data information.               !
!                                                                             !
! Input:                                                                      !
!    NUNIT:  C Unit to read from.                                             !
!    NDATA:  Size of array DATA (Should be >= NDAT as computed herein.)       !
!                                                                             !
! Output:                                                                     !
!    DATA:  The unpacked data array                                           !
!    IERR:  Error flag, non-zero means there was a problem.                   !
!                                                                             !
! Side Effects:                                                               !
!    * Header arrays SEC0, SEC1, SEC2, SEC3, SEC4, XEC4, INFOGRID and         !
!      INFOGRID are filled.                                                   !
!    * The BITMAP array is filled.                                            !
!    * The C unit pointer is advanced to the end of the GRIB record.          !
!                                                                             !
! Modules:                                                                    !
!      MODULE_GRIB                                                            !
!                                                                             !
! Externals:                                                                  !
!      GRIBGET                                                                !
!      GRIBHEADER                                                             !
!      GRIBDATA                                                               !
!                                                                             !
!-----------------------------------------------------------------------------!
  use module_grib

  implicit none

  integer :: nunit
  integer :: debug_level
  integer :: ierr
  real, allocatable, dimension(:) :: datarray
  integer :: ndata
  real, dimension(ndata) :: data

  integer :: ni, nj

  ierr = 0

  call gribget(nunit, ierr)
  if (ierr.ne.0) return

! Unpack the header information

  call gribheader(debug_level,ierr)

! Determine the size of the data array from the information in the header, 
! and allocate the array DATARRAY to hold that data.

  if (sec2(4).ne.50) then
     ni = infogrid(1)
     nj = infogrid(2)
     allocate(datarray(ni*nj))
  else
     ni = (infogrid(1)+1) * (infogrid(1)+2)
     nj = 1
     allocate(datarray(ni*nj))
  endif

! Unpack the data from the GRIB record, and fill the array DATARRAY.

  call gribdata(datarray, ni*nj) 

  data(1:ni*nj) = datarray(1:ni*nj)
#if defined (CRAY)
#else
  deallocate(grec, datarray)
#endif

end subroutine gribread
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine get_sec1(ksec1)
! Return the GRIB Section 1 header information, which has already been
! unpacked by subroutine GRIBHEADER.
  use module_grib
  integer, dimension(100) :: ksec1
  ksec1 = sec1
end subroutine get_sec1
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine get_sec2(ksec2)
! Return the GRIB Section 2 header information, which has already been
! unpacked by subroutine GRIBHEADER.
  use module_grib
  integer, dimension(10) :: ksec2
  ksec2 = sec2
end subroutine get_sec2
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine get_gridinfo(iginfo, ginfo)
  use module_grib
  integer, dimension(40) :: iginfo
  real, dimension(40) :: ginfo
  iginfo = infogrid
  ginfo = gridinfo
end subroutine get_gridinfo
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribprint(isec)
  use module_grib
  implicit none
  integer :: isec
  integer :: ou = 6
  character(len=12) :: string = ',t45,":",i8)'
  character(len=15) :: rstring = ',t45,":",f12.5)'

  if (isec.eq.0) then
     write(*,'(/,"GRIB SECTION 0:")')
     write(ou,'(5x,"Grib Length"'//string) sec0(1)
     write(ou,'(5x,"Grib Edition"'//string) sec0(2)
  else if (isec.eq.1) then
     write(*,'(/,"GRIB SECTION 1:")')
     write(ou,'(5x,"Length of PDS"'//string) sec1(1)
     write(ou,'(5x,"Parameter Table Version"'//string) sec1(2)
     write(ou,'(5x,"Center ID"'//string) sec1(3)
     write(ou,'(5x,"Process ID"'//string) sec1(4)
     write(ou,'(5x,"Grid ID"'//string) sec1(5)
     if (sec1(25) == 1) then
        write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,":     Yes")')
     else if (sec1(25) == 0) then
        write(ou,'(5x,"Is there a Grid Desc. Section (GDS)?",t45,":      No")')
     else
        print*, 'Unrecognized sec1(25): ', sec1(25)
     endif
     if (sec1(26) == 1) then
        write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,":     Yes")')
     else if (sec1(26) == 0) then
        write(ou,'(5x,"Is there a Bit Map Section (BMS)?",t45,":      No")')
     else
        print*, 'Unrecognized sec1(26): ', sec1(26)
     endif
     write(ou,'(5x,"Parameter"'//string) sec1(7)
     write(ou,'(5x,"Level type"'//string) sec1(8)
     if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
          (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
          (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
          (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) ) then
        write(ou,'(5x,"Hgt, pres, etc. of layer top "'//string) sec1(9)
        write(ou,'(5x,"Hgt, pres, etc. of layer bottom "'//string) sec1(10)
     else
        write(ou,'(5x,"Height, pressure, etc "'//string) sec1(9)
     endif
     write(ou,'(5x,"Year"'//string) sec1(11)
     write(ou,'(5x,"Month"'//string) sec1(12)
     write(ou,'(5x,"Day"'//string) sec1(13)
     write(ou,'(5x,"Hour"'//string) sec1(14)
     write(ou,'(5x,"Minute"'//string) sec1(15)
     write(ou,'(5x,"Forecast time unit"'//string) sec1(16)
     write(ou,'(5x,"P1"'//string) sec1(17)
     write(ou,'(5x,"P2"'//string) sec1(18)
     write(ou,'(5x,"Time Range Indicator"'//string) sec1(19)
     write(ou,'(5x,"Number in Ave?"'//string) sec1(20)
     write(ou,'(5x,"Number missing from ave?"'//string) sec1(21)
     write(ou,'(5x,"Century"'//string) sec1(22)
     write(ou,'(5x,"Sub-center"'//string) sec1(23)
     write(ou,'(5x,"Decimal scale factor"'//string) sec1(24)
  elseif ((isec.eq.2) .and. ((sec1(6).eq.128).or.(sec1(6).eq.192))) then
     write(*,'(/,"GRIB SECTION 2:")')
     write(ou,'(5x,"Length of GRID Desc. Section"'//string) sec2(1)
     if ((sec2(2) /= 0).or.(sec2(3) /= 0) .or. (sec2(4) /= 0)) then
        write(ou,'(5x,"Number of V. Coordinate Parms"'//string) sec2(2)
        write(ou,'(5x,"List Starting point"'//string) sec2(3)
        write(ou,'(5x,"Data Representation type"'//string) sec2(4)
     endif

     if (sec2(4).eq.0) then
        write(ou,'(5x,"Cylindrical Equidistant Grid")')
        write(ou,'(10x,"NI"'//string) infogrid(1)
        write(ou,'(10x,"NJ"'//string) infogrid(2)
        write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
        write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
        write(ou,'(10x,"Resolution and Component:", t45,":",B8.8)') infogrid(5)
        write(ou,'(10x,"Lat NI"'//string) infogrid(6)
        write(ou,'(10x,"Lon NJ"'//string) infogrid(7)
        write(ou,'(10x,"Delta-Lon"'//string) infogrid(8)
        write(ou,'(10x,"Delta-Lat"'//string) infogrid(9)
        write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
        write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
        write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)

     else if (sec2(4).eq.1) then
        write(ou,'(5x,"Mercator Grid")')
        write(ou,'(10x,"NI"'//string) infogrid(1)
        write(ou,'(10x,"NJ"'//string) infogrid(2)
        write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
        write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
        write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
        write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
        write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
        write(ou,'(10x,"Dx"'//rstring) gridinfo(8)
        write(ou,'(10x,"Dy"'//rstring) gridinfo(9)
        write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
        write(ou,'(10x,"Latin"'//rstring) gridinfo(11)
        write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
        write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)

     else if (sec2(4).eq.4) then
        write(ou,'(5x,"Gaussian Grid")')
        write(ou,'(10x,"NI"'//string) infogrid(1)
        write(ou,'(10x,"NJ"'//string) infogrid(2)
        write(ou,'(10x,"Original (stored) Lat 1"'//rstring) gridinfo(18)
        write(ou,'(10x,"Lat 1"'//rstring) gridinfo(3)
        write(ou,'(10x,"Lon 1"'//rstring) gridinfo(4)
        write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
        write(ou,'(10x,"Original (stored) Lat NI"'//rstring) gridinfo(17)
        write(ou,'(10x,"Lat NI"'//rstring) gridinfo(6)
        write(ou,'(10x,"Lon NJ"'//rstring) gridinfo(7)
        write(ou,'(10x,"Delta-Lon"'//rstring) gridinfo(8)
        write(ou,'(10x,"Delta-Lat"'//rstring) gridinfo(19)
        write(ou,'(10x,"Number of lats (pole - eq)"'//string) infogrid(9)
        write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
        write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
        write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
     elseif (sec2(4).eq.3) then
        write(ou,'(5x,"Lambert Conformal Grid")')
        write(ou,'(10x,"NI"'//string) infogrid(1)
        write(ou,'(10x,"NJ"'//string) infogrid(2)
        write(ou,'(10x,"Lat 1"'//string) infogrid(3)
        write(ou,'(10x,"Lon 1"'//string) infogrid(4)
        write(ou,'(10x,"Resolution and Component",t45,":", B8.8)') infogrid(5)
        write(ou,'(10x,"Lov"'//string) infogrid(6)
        write(ou,'(10x,"Dx"'//string) infogrid(7)
        write(ou,'(10x,"Dy"'//string) infogrid(8)
        write(ou,'(10x,"Projection center"'//string) infogrid(9)
        write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
        write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
        write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
        write(ou,'(10x,"Latin 1"'//string) infogrid(11)
        write(ou,'(10x,"Latin 2"'//string) infogrid(12)
        write(ou,'(10x,"Lat of southern pole"'//string) infogrid(13)
        write(ou,'(10x,"Lon of southern pole"'//string) infogrid(14)
     elseif (sec2(4).eq.5) then
        write(ou,'(5x,"Polar Stereographic Grid")')
        write(ou,'(10x,"NI"'//string) infogrid(1)
        write(ou,'(10x,"NJ"'//string) infogrid(2)
        write(ou,'(10x,"Lat 1"'//string) infogrid(3)
        write(ou,'(10x,"Lon 1"'//string) infogrid(4)
        write(ou,'(10x,"Resolution and Component", t45,":",B8.8)') infogrid(5)
        write(ou,'(10x,"Lov"'//string) infogrid(6)
        write(ou,'(10x,"Dx"'//string) infogrid(7)
        write(ou,'(10x,"Dy"'//string) infogrid(8)
        write(ou,'(10x,"Projection center"'//string) infogrid(9)
        write(ou,'(10x,"Scanning mode"'//string) infogrid(10)
        write(ou,'(10x,"I-Scanning increment"'//string) infogrid(21)
        write(ou,'(10x,"J-Scanning increment"'//string) infogrid(22)
     elseif (sec2(4).eq.50) then
        write(ou,'(5x,"Spherical harmonic components")')
        write(ou,'(10x,"J-Pentagonal resolution parm:"'//string) infogrid(1)
        write(ou,'(10x,"K-Pentagonal resolution parm:"'//string) infogrid(2)
        write(ou,'(10x,"M-Pentagonal resolution parm:"'//string) infogrid(3)
        write(ou,'(10x,"Representation type"'//string) infogrid(4)
        write(ou,'(10x,"Coefficient storage mode"'//string) infogrid(5)
     endif
  elseif ((isec.eq.3) .and. (sec1(26).eq.1)) then
     write(*,'(/,"GRIB SECTION 3:")')
     write(ou,'(5x,"Length of bit map section"'//string) sec3(1)
     write(ou,'(5x,"Number of unused bits"'//string) sec3(2)
     write(ou,'(5x,"Numeric"'//string) sec3(3)

  elseif (isec.eq.4) then
     write(*,'(/,"GRIB SECTION 4:")')
     write(ou,'(5x,"Length of BDS"'//string) sec4(1)
     write(ou,'(5x,"0/1: grid-point or sph. harm. data"'//string) sec4(2)
     write(ou,'(5x,"0/1: simple or complex packing"'//string) sec4(3)
     write(ou,'(5x,"0/1: floating or integer"'//string) sec4(4)
     write(ou,'(5x,"0/1: No addl flags or addl flags"'//string) sec4(5)
     write(ou,'(5x,"Unused bits"'//string) sec4(6)
     write(ou,'(5x,"Binary Scale Factor"'//string) sec4(7)
     write(ou,'(5x,"Reference Value", t45, ":", F18.8)') xec4(1)
     write(ou,'(5x,"Number of bits for packing"'//string) sec4(8)
  endif

end subroutine gribprint
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine get_bitmap(bm8, ndat)
  use module_grib
  integer, dimension(ndat) :: bm8
  if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
     bm8 = bitmap
  else
     bm8 = 1
  endif
end subroutine get_bitmap
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribxyll(x, y, xlat, xlon)
  use module_grib
  implicit none

  real , intent(in) :: x, y
  real , intent(out) :: xlat, xlon

  real :: r, xkm, ykm, y1
  integer :: iscan, jscan

  if (sec2(4).eq.0) then ! Cylindrical equidistant grid

     xlat = gridinfo(3) + gridinfo(9)*(y-1.)
     xlon = gridinfo(4) + gridinfo(8)*(x-1.)

  elseif (sec2(4) == 1) then      ! Mercator grid
     r = grrth*cosd(gtrue1)
     xkm = (x-1.)*gridinfo(8)
     ykm = (y-1.)*gridinfo(9)
     xlon = gridinfo(4) + (xkm/r)*(180./pi)
     y1 = r*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))/gridinfo(9)
     xlat = 90. - 2. * atan(exp(-gridinfo(9)*(y+y1-1.)/r))*180./pi

  elseif (sec2(4) == 3) then      ! Lambert Conformal grid
     gclon = gridinfo(6)
     r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
     xlat = 90. - 2.*atand(tand(45.-gtrue1/2.)* &
          ((r*gkappa*gridinfo(7))/(grrth*sind(90.-gtrue1)))**(1./gkappa))
     xlon = atan2d((x-1.+gx1),-(y-1.+gy1))/gkappa + gclon

  elseif (sec2(4) == 5) then  ! Polar Stereographic grid
     gclon = gridinfo(6)
     r = sqrt((x-1.+gx1)**2 + (y-1+gy1)**2)
     xlat = 90. - 2.*atan2d((r*gridinfo(7)),(grrth*(1.+sind(gtrue1))))
     xlon = atan2d((x-1.+gx1),-(y-1.+gy1)) + gclon

  elseif (sec2(4) == 4) then ! Gaussian grid

     xlon = gridinfo(4) + gridinfo(8)*(x-1.)
     xlat = gridinfo(3) + gridinfo(19)*(y-1.)

  else
     write(*,'("Unrecognized projection:", I10)') sec2(4)
     write(*,'("STOP in GRIBXYLL")')
     stop
  endif

end subroutine gribxyll
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribllxy(xlat, xlon, x, y)
  use module_grib
  implicit none
  real , intent(in) :: xlat, xlon
  real , intent(out) :: x, y

  real :: r, y1

  if (sec2(4) == 0) then      ! Cylindrical Equidistant grid

     x = 1. + (xlon-gridinfo(4)) / gridinfo(9)
     y = 1. + (xlat-gridinfo(3)) / gridinfo(8)

  else if (sec2(4) == 1) then      ! Mercator grid

     r = grrth*cosd(gtrue1)
     x = 1.+( (r/gridinfo(8)) * (xlon-gridinfo(4)) * (pi/180.) )
     y1 = (r/gridinfo(9))*alog((1.+sind(gridinfo(3)))/cosd(gridinfo(3)))
     y = 1. + ((r/gridinfo(9))*alog((1.+sind(xlat))/cosd(xlat)))-y1

  else if (sec2(4) == 3) then      ! Lambert Conformal grid
     gclon = gridinfo(6)
     r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
          (tand(45.-xlat/2.)/tand(45.-gtrue1/2.)) ** gkappa
     x =  r*sind(gkappa*(xlon-gclon)) - gx1 + 1.
     y = -r*cosd(gkappa*(xlon-gclon)) - gy1 + 1.

  elseif (sec2(4) == 5) then  ! Polar Stereographic grid
     gclon = gridinfo(6)
     r = grrth/gridinfo(7) * tand((90.-xlat)/2.) * (1.+sind(gtrue1))
     x = ( r * sind(xlon-gclon)) - gx1 + 1.
     y = (-r * cosd(xlon-gclon)) - gy1 + 1.

  else
     write(*,'("Unrecognized projection:", I10)') sec2(4)
     write(*,'("STOP in GRIBLLXY")')
     stop
  endif

end subroutine gribllxy
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine glccone (fsplat,ssplat,sign,confac)
  use module_grib
  implicit none
  real, intent(in) :: fsplat,ssplat
  integer, intent(in) :: sign
  real, intent(out) :: confac
  if (abs(fsplat-ssplat).lt.1.E-3) then
     confac = sind(fsplat)
  else
     confac = log10(cosd(fsplat))-log10(cosd(ssplat))
     confac = confac/(log10(tand(45.-float(sign)*fsplat/2.))- &
          log10(tand(45.-float(sign)*ssplat/2.)))
  endif
end subroutine glccone
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
subroutine gribheader(debug_level,ierr)
!
! IERR non-zero means there was a problem unpacking the grib header.
!
  use module_grib
  implicit none
  integer :: debug_level
  integer :: ierr

  integer, parameter :: nsec1 = 24

  integer, dimension(nsec1) :: &
       iw1=(/3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,2/)
  integer :: icount, iskip, ibts, nbm, nbz, i9skip, i17skip

  integer :: iman, ichar, isign, iscan

  integer,  allocatable, dimension(:) :: bm8

  real :: r
  integer :: isvns
  integer :: gsvns = 0

  if (gsvns.eq.0) then
     if (mwsize.eq.32) then
        gsvns = transfer('7777', gsvns)
     elseif(mwsize.eq.64) then
        call gbyte_g1(char(0)//char(0)//char(0)//char(0)//'7777', gsvns, 0, mwsize)
     endif
  endif

! Section 0:
  sec0(2) = ied
  if (ied.eq.1) then
     call gbyte_g1(grec, sec0(1), 32, 24)
     iskip = 64
  elseif (ied.eq.0) then
     sec0(1) = gribsize(grec,200, ierr)
     iskip = 32
  endif

! Section 1:
  i9skip = iskip + 80
  i17skip = iskip + 144
  do icount = 1, nsec1 - ((1-ied)*6)
     ibts = iw1(icount)*8
     call gbyte_g1(grec, sec1(icount), iskip, ibts)
     iskip = iskip + ibts
  enddo
  if (ied.eq.0) sec1(22) = 20
  ! Sec1 indices 9 and 10 might actually be one value, not two.
  ! If this is the case, reread sec1(9), and set sec1(10) to zero:
  if ( (sec1(8) == 101) .or. (sec1(8) == 104) .or. (sec1(8) == 106) .or. &
       (sec1(8) == 108) .or. (sec1(8) == 110) .or. (sec1(8) == 112) .or. &
       (sec1(8) == 114) .or. (sec1(8) == 116) .or. (sec1(8) == 120) .or. &
       (sec1(8) == 121) .or. (sec1(8) == 128) .or. (sec1(8) == 141) .or. &
       (sec1(8) == 236) ) then
     ! No action here.
  else
     call gbyte_g1(grec, sec1(9), i9skip, 16)
     sec1(10) = 0.
  endif

  if (sec1(24).ge.32768) sec1(24) = 32768-sec1(24)

  ! If the TIME/RANGE INDICATOR (sec1(19)) indicates that the time P1
  ! is spread over two bytes, then recompute sec1(17) and set sec1(18)
  ! to zero.
  if (sec1(19) == 10) then
     call gbyte_g1(grec, sec1(17), i17skip, 16)
     sec1(18) = 0
  endif

  ! Pull out single bits from sec1(6) for the GDS and BMS flags:
  sec1(25) = sec1(6)/128
  sec1(26) = mod(sec1(6)/64,2)

! Section 2:
! if ((sec1(6) == 128) .or. (sec1(6) == 192)) then
  if (sec1(25) == 1) then

     if (ied.eq.0) then
        iskip = 32 + sec1(1)*8
     elseif (ied.eq.1) then
        iskip = 64 + sec1(1)*8
     endif
     call gbyte_g1(grec, sec2(1), iskip, 24)
     iskip = iskip + 24
     call gbytes_g1(grec, sec2(2), iskip, 8, 0, 3)
     iskip = iskip + 8*3

     if (sec2(4) == 0) then
        ! Lat/Lon Grid:
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)        
        iskip = iskip + 32
        call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
        iskip = iskip + 48
        call gbyte_g1(grec, infogrid(5), iskip, 8)
        iskip = iskip + 8
        call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
        iskip = iskip + 48
        call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)        
        iskip = iskip + 32
        call gbyte_g1(grec, infogrid(21), iskip, 1)
        infogrid(21) = 1-(infogrid(21)*2)
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(22), iskip, 1)
        infogrid(22) = (infogrid(22)*2)-1
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(10), iskip, 1)
        iskip = iskip + 1
        iskip = iskip + 5
        call gbyte_g1(grec, infogrid(11), iskip, 32)
        iskip = iskip + 32

!MGD        if ( debug_level .gt. 100 ) then
!MGD        print *, "lat/lon grib grid info", infogrid(1), infogrid(3), &
!MGD        infogrid(5), infogrid(6), infogrid(8), infogrid(21), &
!MGD        infogrid(22), infogrid(10), infogrid(11), infogrid(8)
!MGD        end if

        infogrid(8) = infogrid(8) * infogrid(21)
        infogrid(9) = infogrid(9) * infogrid(22)

        gridinfo(1) = float(infogrid(1))
        gridinfo(2) = float(infogrid(2))
        if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
        if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
        gridinfo(3) = float(infogrid(3))*0.001
        gridinfo(4) = infogrid(4) * 0.001
        if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
        if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
        gridinfo(6) = infogrid(6) * 0.001
        gridinfo(7) = infogrid(7) * 0.001
        gridinfo(8) = infogrid(8) * 0.001
        gridinfo(9) = infogrid(9) * 0.001
        gridinfo(21) = float(infogrid(21))
        gridinfo(22) = float(infogrid(22))
     elseif (sec2(4) == 1) then ! Mercator grid
        ! Number of points in X and Y
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)
        iskip = iskip + 32
        ! Starting lat and lon
        call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
        iskip = iskip + 48
        ! "Resolution and component flags"
        call gbyte_g1(grec, infogrid(5), iskip, 8)
        iskip = iskip + 8
        ! Ending lat and lon
        call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
        iskip = iskip + 48
        ! Truelat, 3 bytes
        call gbyte_g1(grec, infogrid(11), iskip, 24)
        iskip = iskip + 24
        ! "Reserved", i.e., skip a byte
        iskip = iskip + 8
        ! Scanning mode flags, first three bits of the next byte
        ! and skip the last five bits.
        call gbyte_g1(grec, infogrid(21), iskip, 1)
        infogrid(21) = 1-(infogrid(21)*2)
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(22), iskip, 1)
        infogrid(22) = (infogrid(22)*2)-1
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(10), iskip, 1)
        iskip = iskip + 1
        iskip = iskip + 5
        ! Grid increment in X and Y
        call gbytes_g1(grec, infogrid(8), iskip, 24, 0, 2)
        iskip = iskip + 48
        ! Done reading map specifications.
        ! Now do various conversions:

        gridinfo(1) = float(infogrid(1)) ! ok
        gridinfo(2) = float(infogrid(2)) ! ok

        if (infogrid(3) .ge.8388608) infogrid(3)  = 8388608 - infogrid(3)
        if (infogrid(4) .ge.8388608) infogrid(4)  = 8388608 - infogrid(4)
        if (infogrid(6) .ge.8388608) infogrid(6)  = 8388608 - infogrid(6)
        if (infogrid(7) .ge.8388608) infogrid(7)  = 8388608 - infogrid(7)
        if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
        gridinfo(3)  = infogrid(3)  * 0.001
        gridinfo(4)  = infogrid(4)  * 0.001
        gridinfo(6)  = infogrid(6)  * 0.001
        gridinfo(7)  = infogrid(7)  * 0.001
        gridinfo(8)  = infogrid(8)  * 0.001
        gridinfo(9)  = infogrid(9)  * 0.001
        gridinfo(11) = infogrid(11) * 0.001

        gridinfo(21) = infogrid(21)
        gridinfo(22) = infogrid(22)

        gridinfo(20) = 6370.949
        grrth = gridinfo(20)
        gtrue1 = gridinfo(11)

     elseif (sec2(4) == 3) then
        if (ied.eq.0) then
           print '(//,"*** Despair ***"//)'
           stop
        endif
! Lambert Conformal:
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)        
        iskip = iskip + 32
        call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
        iskip = iskip + 48
        if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
        if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
        call gbyte_g1(grec, infogrid(5), iskip, 8)
        iskip = iskip + 8
        call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3)
        if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
        iskip = iskip + 72
        call gbyte_g1(grec, infogrid(9), iskip, 8)
        iskip = iskip + 8
        call gbyte_g1(grec, infogrid(21), iskip, 1)
        infogrid(21) = 1-(infogrid(21)*2)
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(22), iskip, 1)
        infogrid(22) = (infogrid(22)*2)-1
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(10), iskip, 1)
        iskip = iskip + 1
        iskip = iskip + 5
        call gbytes_g1(grec, infogrid(11), iskip, 24, 0, 4)
        if (infogrid(11).ge.8388608) infogrid(11) = 8388608 - infogrid(11)
        if (infogrid(12).ge.8388608) infogrid(12) = 8388608 - infogrid(12)
        if (infogrid(13).ge.8388608) infogrid(13) = 8388608 - infogrid(13)
        if (infogrid(14).ge.8388608) infogrid(14) = 8388608 - infogrid(14)
        iskip = iskip + 96
        call gbyte_g1(grec, infogrid(15), iskip, 16)
        iskip = iskip + 16

        infogrid(7) = infogrid(7) * infogrid(21)
        infogrid(8) = infogrid(8) * infogrid(22)


        gridinfo(1) = float(infogrid(1))
        gridinfo(2) = float(infogrid(2))
        gridinfo(3) = infogrid(3) * 0.001
        gridinfo(4) = infogrid(4) * 0.001
        gridinfo(6) = infogrid(6) * 0.001
        gridinfo(7) = infogrid(7) * 0.001
        gridinfo(8) = infogrid(8) * 0.001
        gridinfo(9) = infogrid(9) * 0.001
        gridinfo(11) = infogrid(11) * 0.001
        gridinfo(12) = infogrid(12) * 0.001
        gridinfo(13) = infogrid(13) * 0.001
        gridinfo(14) = infogrid(14) * 0.001

        gridinfo(20) = 6370
        ! a priori knowledge:
        if (sec1(5).eq.212) then
           gridinfo(3) = 12.190
           gridinfo(4) = -133.459
           gridinfo(7) = 40.63525
           gridinfo(8) = 40.63525
           gridinfo(20) = 6370
        endif

!=============================================================================!
! More a priori knowledge:                                                    !
! Correct some bad lat/lon numbers coded into some RUC headers.               !
!                                                                             !
        if (sec1(3) == 59) then  ! If FSL
           if (sec1(4) == 86) then  ! and RUC
              if (sec1(5) == 255) then 
           ! Check to correct bad lat/lon numbers.
                 if (infogrid(3) == 16909) then
                    infogrid(3) = 16281
                    gridinfo(3) = 16.281
                 endif
                 if (infogrid(4) == 236809) then
                    infogrid(4) = 2338622
                    gridinfo(4) = 233.8622
                 endif
              endif
           endif
        endif
!=============================================================================!


        gridinfo(21) = float(infogrid(21))
        gridinfo(22) = float(infogrid(22))

        ! Map parameters
        glat1 = gridinfo(3)
        glon1 = gridinfo(4)
        gclon = gridinfo(6)
        if (gclon.gt.180.) gclon = -(360.-gclon)
        if ((gclon<0).and.(glon1>180)) glon1 = glon1-360.  
        gtrue1 = gridinfo(11)
        gtrue2 = gridinfo(12)
        grrth = gridinfo(20)
        call glccone(gtrue1, gtrue2, 1, gkappa)
        r = grrth/(gridinfo(7)*gkappa)*sind(90.-gtrue1) * &
             (tand(45.-glat1/2.)/tand(45.-gtrue1/2.)) ** gkappa
        gx1 =  r*sind(gkappa*(glon1-gclon))
        gy1 = -r*cosd(gkappa*(glon1-gclon))

     elseif (sec2(4) == 4) then
        ! Gaussian Grid:
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)        
        iskip = iskip + 32
        call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)
        iskip = iskip + 48
        call gbyte_g1(grec, infogrid(5), iskip, 8)
        iskip = iskip + 8
        call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 2)
        iskip = iskip + 48
        call gbytes_g1(grec, infogrid(8), iskip, 16, 0, 2)        
        iskip = iskip + 32
        call gbyte_g1(grec, infogrid(21), iskip, 1)
        infogrid(21) = 1-(infogrid(21)*2)
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(22), iskip, 1)
        infogrid(22) = (infogrid(22)*2)-1
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(10), iskip, 1)
        iskip = iskip + 1
        iskip = iskip + 5
        call gbyte_g1(grec, infogrid(11), iskip, 32)
        iskip = iskip + 32

        infogrid(8) = infogrid(8) * infogrid(21)

        gridinfo(1) = float(infogrid(1))
        gridinfo(2) = float(infogrid(2))
        if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
        if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
        gridinfo(3) = float(infogrid(3))*0.001
        gridinfo(4) = infogrid(4) * 0.001
        if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)
        if (infogrid(7).ge.8388608) infogrid(7) = 8388608 - infogrid(7)
        gridinfo(6) = infogrid(6) * 0.001
        gridinfo(7) = infogrid(7) * 0.001
        gridinfo(8) = infogrid(8) * 0.001
        gridinfo(21) = float(infogrid(21))
        gridinfo(22) = float(infogrid(22))

        ! Compute an approximate delta-latitude and starting latitude.
        ! Replace the stored value of starting latitude with approximate one.
        gridinfo(18) = gridinfo(3)
        infogrid(18) = infogrid(3)
        gridinfo(17) = gridinfo(6)
        infogrid(17) = infogrid(6)
!        call griblgg(infogrid(2), gridinfo(3), gridinfo(19))
!        infogrid(19) = nint(gridinfo(19)*1000.)
!        infogrid(3) = nint(gridinfo(3)*1000.)
        gridinfo(6) = -gridinfo(3)
        infogrid(6) = -infogrid(3)

     elseif (sec2(4) == 5) then
! Polar Stereographic Grid
        if (ied.eq.0) then
           print '(//,"*** Despair ***"//)'
           stop
        endif
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 2)  ! NX and NY
        iskip = iskip + 32
        call gbytes_g1(grec, infogrid(3), iskip, 24, 0, 2)  ! LAT1 and LON1
        iskip = iskip + 48
        call gbyte_g1(grec, infogrid(5), iskip, 8) ! Resolution and Component
        iskip = iskip + 8
        call gbytes_g1(grec, infogrid(6), iskip, 24, 0, 3) ! LOV, DX, and DY
        iskip = iskip + 72
        call gbyte_g1(grec, infogrid(9), iskip, 8) ! Projection center flag
        iskip = iskip + 8
        call gbyte_g1(grec, infogrid(21), iskip, 1)
        infogrid(21) = 1-(infogrid(21)*2)
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(22), iskip, 1)
        infogrid(22) = (infogrid(22)*2)-1
        iskip = iskip + 1
        call gbyte_g1(grec, infogrid(10), iskip, 1)
        iskip = iskip + 1
        iskip = iskip + 5
!         call gbyte_g1(grec, infogrid(11), iskip, 32) ! Set to 0 (reserved)
        iskip = iskip + 32

        if (infogrid(3).ge.8388608) infogrid(3) = 8388608 - infogrid(3)
        if (infogrid(4).ge.8388608) infogrid(4) = 8388608 - infogrid(4)
        if (infogrid(6).ge.8388608) infogrid(6) = 8388608 - infogrid(6)


        infogrid(7) = infogrid(7) * infogrid(21)
        infogrid(8) = infogrid(8) * infogrid(22)

        gridinfo(1) = float(infogrid(1))
        gridinfo(2) = float(infogrid(2))
        gridinfo(3) = infogrid(3) * 0.001
        gridinfo(4) = infogrid(4) * 0.001
        gridinfo(6) = infogrid(6) * 0.001
        gridinfo(7) = infogrid(7) * 0.001
        gridinfo(8) = infogrid(8) * 0.001

        gridinfo(20) = 6370

        ! a priori knowledge:
        if (sec1(5).eq.240) then
           gridinfo(3) = 22.7736
           gridinfo(4) = -120.376
           gridinfo(7) = 4.7625
           gridinfo(8) = 4.7625
           gridinfo(20) = 6370
        endif

        ! Map parameters
        glat1 = gridinfo(3)
        glon1 = gridinfo(4)
        gclon = gridinfo(6)
        if (gclon.gt.180.) gclon = -(360.-gclon)
        ! GRIB edition 1 Polar Stereographic grids are true at 60 degrees
        ! Which hemisphere depends on infogrid(9), the "Projection Center Flag"
        grrth = gridinfo(20)
        if (infogrid(9) > 127) then
           gtrue1 = -60. 
           r = grrth/gridinfo(7) * tand((-90.-glat1)/2.) * (1.+sind(-gtrue1))
           gx1 = -r * sind(glon1-gridinfo(6))
           gy1 = -r * cosd(glon1-gridinfo(6)) 
        else
           gtrue1 = 60. 
           r = grrth/gridinfo(7) * tand((90.-glat1)/2.) * (1.+sind(gtrue1))
           gx1 = r * sind(glon1-gridinfo(6))
           gy1 = -r * cosd(glon1-gridinfo(6))
        endif

        gridinfo(21) = float(infogrid(21))
        gridinfo(22) = float(infogrid(22))

     elseif (sec2(4) == 50) then
! Spherical harmonic coefficients
        if (ied.eq.0) then
           print '(//,"*** Despair ***"//)'
           stop
        endif
        call gbytes_g1(grec, infogrid(1), iskip, 16, 0, 3)
        iskip = iskip + 48
        call gbytes_g1(grec, infogrid(4), iskip, 8, 0, 2)
        iskip = iskip + 16

        iskip = iskip + 18*8

     else
        call gribprint(0)
        call gribprint(1)
        call gribprint(2)
        call gribprint(3)
        call gribprint(4)
        write(*,'("Unrecognized grid: ", i8)') sec2(4)
        write(*,'("This grid is not currently supported.")')
        write(*,'("Write your own program to put the data to the intermediate format")')
        stop
     endif

  endif

! Section 3
  if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then
     if (ied.eq.0) then
        print '(//,"*** Despair ***"//)'
        stop
     endif

     if (ied.eq.0) then
        iskip = 32 + sec1(1)*8 + sec2(1)*8
     elseif (ied.eq.1) then
        iskip = 64 + sec1(1)*8 + sec2(1)*8
     endif
     call gbyte_g1(grec, sec3(1), iskip, 24)
     iskip = iskip + 24
     call gbyte_g1(grec, sec3(2), iskip, 8)
     iskip = iskip + 8
     call gbyte_g1(grec, sec3(3), iskip, 16)
     iskip = iskip + 16

#if defined (CRAY)
#else
     allocate(bitmap((sec3(1)-6)*8))
#endif
     allocate(bm8((sec3(1)-6)*8))
     call gbytes_g1(grec, bm8, iskip, 1, 0, (sec3(1)-6)*8)
     bitmap(1:size(bm8)) = bm8(1:size(bm8))
     deallocate(bm8)
     iskip = iskip + sec3(1)-6
  else
     sec3 = 0
  endif

! Section 4
  if ((sec1(6).eq.128).or.(sec1(6).eq.192)) then
     if (ied.eq.0) then
        iskip = 32 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
     elseif (ied.eq.1) then
        iskip = 64 + sec1(1)*8 + sec2(1)*8 + sec3(1)*8
     endif
     call gbyte_g1(grec, sec4(1), iskip, 24)
     if (sec4(1) > (sec0(1) - sec1(1) - sec2(1) - sec3(1) - 4)) then
        write(*,'(/,"*** I have good reason to believe that this GRIB record is")')
        write(*,'("*** corrupted or miscoded.",/)')
        ierr = 1
        return
     endif
     iskip = iskip + 24
     call gbytes_g1(grec, sec4(2), iskip, 1,0,4)
     iskip = iskip + 4
     call gbyte_g1(grec, sec4(6), iskip, 4)
     iskip = iskip + 4
! Get the binary scale factor
     call gbyte_g1(grec, isign, iskip, 1)
     iskip = iskip + 1
     call gbyte_g1(grec, sec4(7), iskip, 15)
     iskip = iskip + 15
     sec4(7) = sec4(7) * (-2*isign+1)
! Get the reference value:
     call gbyte_g1(grec, isign, iskip, 1)
     iskip = iskip + 1
     isign = -2*isign+1
     call gbyte_g1(grec, ichar, iskip, 7)
     iskip = iskip + 7
     call gbyte_g1(grec, iman, iskip, 24)
     iskip = iskip + 24
     if ( iman .ne. 0 ) then
       xec4(1) = float(isign) * (2.**(-24)) * float(iman) *  &
          (16.**(ichar-64))
     else
       xec4(1) = 0.
     endif

     call gbyte_g1(grec,sec4(8), iskip, 8)
! sec4(8) is the number of bits used per datum value.
! If sec4(8) = 255, assume they mean sec4(8) = 0
     if (sec4(8) == 255) sec4(8) = 0
     iskip = iskip + 8
  endif

! Section 5
  call gbyte_g1(grec, isvns, ((sec0(1)-4)*8), 32)
  if (isvns.ne.gsvns) then
     write(*, '("End-of-record mark (7777) not found", 2I10)') isvns
     write(*, '("Sec0(1) = ", I8, i2)') sec0(1), sevencount
     sevencount = sevencount + 1
     if (sevencount > 10) then
        write(*,'(//," *** Found more than 10 consecutive bad GRIB records")')
        write(*,'(" *** Let''s just stop now.",//)')
        write(*,'(" Perhaps the analysis file should have been converted",/,&
             &" from COS-Blocked format?",//)')
        stop
     endif
  else
     sevencount = 0
  endif

  ierr = 0

end subroutine gribheader
!
!=============================================================================!
!=============================================================================!
!=============================================================================!
!
  subroutine gribdata(datarray, ndat) 

!-----------------------------------------------------------------------------!
!                                                                             !
! Read and unpack the data from a GRIB record.                                !
!                                                                             !
! Input:                                                                      !
!    NDAT:  The size of the data array we expect to unpack.                   !
!                                                                             !
! Output:                                                                     !
!    DATARRAY:  The unpacked data from the GRIB record                        !
!                                                                             !
! Side Effects:                                                               !
!    STOP if it cannot unpack the data.                                       !
!                                                                             !
! Externals:                                                                  !
!     SGUP_BITMAP                                                             !
!     SGUP_NOBITMAP                                                           !
!     CSHUP                                                                   !
!                                                                             !
! Modules:                                                                    !
!     MODULE_GRIB                                                             !
!                                                                             !
!-----------------------------------------------------------------------------!
    use module_grib

    implicit none

    integer :: ndat
    real , dimension(ndat) :: datarray
    integer, dimension(ndat) :: IX

    integer :: iskip, nbm

    if (sec4(2) == 0) then ! Grid-point data
       if (sec4(3).eq.0) then ! Simple unpacking
          if ((sec1(6).eq.64).or.(sec1(6).eq.192)) then ! There is a bitmap
             call SGUP_BITMAP(datarray, ndat)
          else
             call SGUP_NOBITMAP(datarray, ndat)
          endif
       else
          write(*,'(//,"***** No complex unpacking of gridpoint data.")')
          write(*,'("***** Option not yet available.",//)')
!         write(*,'("***** Complain to mesouser@ucar.edu",//)')
          stop
       endif
    else
       if (sec4(3).eq.0) then ! Simple unpacking
          write(*,'(//,"***** No simple unpacking of spherical-harmonic coefficients.")')
          write(*,'("***** Option not yet available.",//)')
!         write(*,'("***** Complain to mesouser@ucar.edu",//)')
          stop
       elseif (sec4(3).eq.1) then
          call CSHUP(datarray, ndat)
       endif
    endif

end subroutine gribdata

subroutine deallogrib
! Deallocates a couple of arrays that may be allocated.
  use module_grib
#if defined (CRAY)
#else
  if (allocated(grec)) deallocate(grec)
  if (allocated(bitmap)) deallocate(bitmap)
#endif
end subroutine deallogrib

SUBROUTINE gribLGG( NLAT, startlat, deltalat )


  implicit none
!
!  LGGAUS finds the Gaussian latitudes by finding the roots of the
!  ordinary Legendre polynomial of degree NLAT using Newtons
!  iteration method.
!
!  On entry:
  integer NLAT ! the number of latitudes (degree of the polynomial)
!
!  On exit: for each Gaussian latitude

  double precision, dimension(NLAT) :: LATG ! Latitude

! Approximations to a regular latitude grid:
  real :: deltalat
  real :: startlat

!-----------------------------------------------------------------------

  integer :: iskip = 15
  double precision :: sum1 = 0.
  double precision :: sum2 = 0.
  double precision :: sum3 = 0.
  double precision :: sum4 = 0.
  double precision :: xn

  integer, save :: SAVE_NLAT = -99
  real, save :: save_deltalat = -99.
  real, save :: save_startlat = -99.

  double precision, dimension(nlat) ::  COSC, SINC
  double precision, parameter :: PI = 3.141592653589793
!
!    -convergence criterion for iteration of cos latitude
  double precision, parameter :: XLIM  = 1.0E-14
  integer :: nzero, i, j
  double precision :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d

  if (nlat == save_nlat) then
     deltalat = save_deltalat
     startlat = save_startlat
     return
  endif
!
!    -the number of zeros between pole and equator
  NZERO = NLAT/2
!
!    -set first guess for cos(colat)
  DO I=1,NZERO
     COSC(I) = SIN( (I-0.5)*PI/NLAT + PI*0.5 )
  ENDDO
!
!    -constants for determining the derivative of the polynomial
  FI  = NLAT
  FI1 = FI+1.0
  A   = FI*FI1 / SQRT(4.0*FI1*FI1-1.0)
  B   = FI1*FI / SQRT(4.0*FI*FI-1.0)
!
!    -loop over latitudes, iterating the search for each root
  DO I=1,NZERO
     J=0
!
!       -determine the value of the ordinary Legendre polynomial for
!       -the current guess root
     LOOP30 : DO 
        CALL LGORD( G, COSC(I), NLAT )
!
!       -determine the derivative of the polynomial at this point
        CALL LGORD( GM, COSC(I), NLAT-1 )
        CALL LGORD( GP, COSC(I), NLAT+1 )
        GT = (COSC(I)*COSC(I)-1.0) / (A*GP-B*GM)
!
!       -update the estimate of the root
        DELTA   = G*GT
        COSC(I) = COSC(I) - DELTA
!
!       -if convergence criterion has not been met, keep trying
        J = J+1
        IF( ABS(DELTA).LE.XLIM ) EXIT LOOP30
     ENDDO LOOP30
  ENDDO
!
!  Determine the sin(colat)
  SINC(1:NZERO) = SIN(ACOS(COSC(1:NZERO)))
!
!    -if NLAT is odd, set values at the equator
  IF( MOD(NLAT,2) .NE. 0 ) THEN
     I       = NZERO+1
     SINC(I) = 1.0
     latg(i) = 0.
  END IF

! Set the latitudes.

  latg(1:NZERO) = dacos(sinc(1:NZERO)) * 180. / pi

! Determine the southern hemisphere values by symmetry
  do i = 1, nzero
     latg(nlat-nzero+i) = -latg(nzero+1-i)
  enddo


! Now that we have the true values, find some approximate values.

  xn = float(nlat-iskip*2)
  do i = iskip+1, nlat-iskip
     sum1 = sum1 + latg(i)*float(i)
     sum2 = sum2 + float(i)
     sum3 = sum3 + latg(i)
     sum4 = sum4 + float(i)**2
  enddo

  b = (xn*sum1 - sum2*sum3) / (xn*sum4 - sum2**2)
  a = (sum3 - b * sum2) / xn

  deltalat = sngl(b)
  startlat = sngl(a + b)

  save_nlat = nlat
  save_deltalat = deltalat
  save_startlat = startlat

contains
  SUBROUTINE LGORD( F, COSC, N )
    implicit none
!
!  LGORD calculates the value of an ordinary Legendre polynomial at a
!  latitude.
!
!  On entry:
!     COSC - cos(colatitude)
!     N      - the degree of the polynomial
!
!  On exit:
!     F      - the value of the Legendre polynomial of degree N at
!              latitude asin(COSC)
    double precision :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang
    integer :: n, k

!------------------------------------------------------------------------

    colat = acos(cosc)
    c1 = sqrt(2.0)
    do k=1,n
       c1 = c1 * sqrt( 1.0 - 1.0/(4*k*k) )
    enddo
    fn = n
    ang= fn * colat
    s1 = 0.0
    c4 = 1.0
    a  =-1.0
    b  = 0.0
    do k=0,n,2
       if (k.eq.n) c4 = 0.5 * c4
       s1 = s1 + c4 * cos(ang)
       a  = a + 2.0
       b  = b + 1.0
       fk = k
       ang= colat * (fn-fk-2.0)
       c4 = ( a * (fn-b+1.0) / ( b * (fn+fn-a) ) ) * c4
    enddo
    f = s1 * c1
  end subroutine lgord

END SUBROUTINE GRIBLGG

SUBROUTINE REORDER_IT (a, nx, ny, dx, dy, iorder)

      use module_debug

      implicit none
      integer :: nx, ny, iorder
      integer :: i, j, k, m
      real :: dx, dy
      real, dimension(nx*ny) :: a, z

      if (iorder .eq. 0 .and. dx .gt. 0. .and. dy .lt. 0) return
      k = 0
      call mprintf(.true.,DEBUG, &
        "Reordering GRIB array : dx = %f  , dy = %f  , iorder = %i",  &
	 f1=dx,f2=dy,i1=iorder)
      if (iorder .eq. 0 ) then
	if ( dx .lt. 0 .and. dy .lt. 0. ) then
	  do j = 1, ny
	  do i = nx, 1, -1
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        else if ( dx .lt. 0 .and. dy .gt. 0. ) then
	  do j = ny, 1, -1
	  do i = nx, 1, -1
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        else if ( dx .gt. 0 .and. dy .gt. 0. ) then
	  do j = ny, 1, -1
	  do i = 1, nx
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        endif
      else
	if ( dx .gt. 0 .and. dy .lt. 0. ) then
	  do i = 1, nx
	  do j = 1, ny
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        else if ( dx .lt. 0 .and. dy .lt. 0. ) then
	  do i = nx, 1, -1
	  do j = 1, ny
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        else if ( dx .lt. 0 .and. dy .lt. 0. ) then
	  do i = nx, 1, -1
	  do j = ny, 1, -1
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        else if ( dx .gt. 0 .and. dy .gt. 0. ) then
	  do i = 1, nx
	  do j = ny, 1, -1
	    k = k + 1
	    m = i * j
	    z(k) = a(m)
	  enddo
	  enddo
        endif
      endif
!  now put it back in the 1-d array and reset the dx and dy
      do k = 1, nx*ny
        a(k) = z(k)
      enddo
      dx = abs ( dx)
      dy = -1 * abs(dy)
      return
END SUBROUTINE REORDER_IT