/WPS/ungrib/src/gribcode.F

Large files files are truncated, but you can click here to view the full file

!                                                                             !
!*****************************************************************************!
!                                                                             !
! 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) inf…