!> @file
!! @brief Determine whether file is grib or not.
!! @author gayno org: w/np2 @date 2007-nov-28

!> Determine whether file is grib or not.
!!
!! program history log:
!! - 2007-nov-28  gayno    - initial version
!! - 2011-apr-26  gayno    - replace my simple-minded logic
!!                         with call to w3lib routin skgb.
!! - 2014-feb-07  gayno    - determine whether file is
!!                         grib1 or grib2.
!!
!! @param[in] file_name - file to be checked
!! @param[out] isgrib - '1' or '2' if grib1/2 file '0' if not grib
!!
!! input files: 
!!     - file to be checked, fort.11
!!
!! condition codes:  all fatal
!!     - bad file open, fort.11
!!          
 subroutine grib_check(file_name, isgrib)

 implicit none

 character*(*), intent(in)         :: file_name
 integer, parameter                :: iunit=11
 integer                           :: istat, iseek, mseek, lskip, lgrib, version
 integer, intent(out)              :: isgrib

 print*,"- CHECK FILE TYPE OF: ", trim(file_name)
 call baopenr (iunit, file_name, istat)

 if (istat /= 0) then
   print*,'- FATAL ERROR: BAD FILE OPEN. ISTAT IS ',istat
   call w3tage('SNOW2MDL')
   call errexit(40)
 end if
 
 iseek = 0
 mseek = 64
 call skgb2(iunit, iseek, mseek, lskip, lgrib, version)

 call baclose(iunit, istat)
 
 if (lgrib > 0) then
   isgrib = version
   if (isgrib == 1) print*,"- FILE IS GRIB1"
   if (isgrib == 2) print*,"- FILE IS GRIB2"
 else
   isgrib = 0
   print*,"- FILE IS BINARY"
 endif

 return

 end subroutine grib_check

!> Determine whether file is grib or not.
!!
!! Based on w3nco library routine skgb.
!!  
!! @param[in] lugb file unit number
!! @param[in] iseek number of bits to skip before search.
!! @param[in] mseek max number of bytes to search.
!! @param[out] lskip number of bytes to skip before message
!! @param[out] lgrib number of bytes in message. '0' if not grib.
!! @param[out] i1 '1' or '2' if grib1/2 file. '0' if not grib.
!!
!! input file:
!! - file to be checked, unit=lugb
!!
!! @author George Gayno org: w/np2 @date 2014-Feb-07
 SUBROUTINE SKGB2(LUGB,ISEEK,MSEEK,LSKIP,LGRIB,I1)
 implicit none
 INTEGER, INTENT( IN)     :: LUGB, ISEEK, MSEEK
 INTEGER, INTENT(OUT)     :: LSKIP, LGRIB, I1
 INTEGER, PARAMETER       :: LSEEK=128
 INTEGER                  :: K, KZ, KS, KG, KN, KM, I4, K4
 CHARACTER Z(LSEEK)
 CHARACTER Z4(4)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 I1=0
 LGRIB=0
 KS=ISEEK
 KN=MIN(LSEEK,MSEEK)
 KZ=LSEEK
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  LOOP UNTIL GRIB MESSAGE IS FOUND
 DO WHILE(LGRIB.EQ.0.AND.KN.GE.8.AND.KZ.EQ.LSEEK)
!  READ PARTIAL SECTION
   CALL BAREAD(LUGB,KS,KN,KZ,Z)
   KM=KZ-8+1
   K=0
!  LOOK FOR 'GRIB...1' IN PARTIAL SECTION
   DO WHILE(LGRIB.EQ.0.AND.K.LT.KM)
     CALL GBYTEC(Z,I4,(K+0)*8,4*8)
     CALL GBYTEC(Z,I1,(K+7)*8,1*8)
     IF(I4.EQ.1196575042.AND.(I1.EQ.1.OR.I1.EQ.2)) THEN
!  LOOK FOR '7777' AT END OF GRIB MESSAGE
       IF (I1.EQ.1) CALL GBYTEC(Z,KG,(K+4)*8,3*8)
       IF (I1.EQ.2) CALL GBYTEC(Z,KG,(K+12)*8,4*8)
       CALL BAREAD(LUGB,KS+K+KG-4,4,K4,Z4)
       IF(K4.EQ.4) THEN
         CALL GBYTEC(Z4,I4,0,4*8)
         IF(I4.EQ.926365495) THEN
!  GRIB MESSAGE FOUND
           LSKIP=KS+K
           LGRIB=KG
         ENDIF
       ENDIF
     ENDIF
     K=K+1
   ENDDO
   KS=KS+KM
   KN=MIN(LSEEK,ISEEK+MSEEK-KS)
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 RETURN
 END subroutine skgb2

!> Convert from the grib2 grid description template array
!! used by the ncep grib2 library, to the grib1 grid
!! description section array used by ncep ipolates library.
!!
!! @param[in] igdtnum grib2 grid desc template number
!! @param[in] igdstmpl grib2 grid desc template array
!! @param[in] igdtlen grib2 grid desc template array size
!! @param[out] kgds grib1 grid description section array used by ncep ipolates library.
!! @param[out] ni i grid dimensions
!! @param[out] nj j grid dimensions
!! @param[out] res grid resolution in km
!!
!! condition codes:
!!   50 - unrecognized model grid type; fatal
!!
!! @author George Gayno org: w/np2 @date 2014-Sep-26
 subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)

 implicit none

 integer, intent(in   )  :: igdtnum, igdtlen, igdstmpl(igdtlen)
 integer, intent(  out)  :: kgds(200), ni, nj
 integer                 :: iscale

 real,    intent(  out)  :: res

 kgds=0

 if (igdtnum.eq.0) then        ! lat/lon grid

   iscale=igdstmpl(10)*igdstmpl(11)
   if (iscale == 0) iscale = 1e6
   kgds(1)=0                   ! oct 6
   kgds(2)=igdstmpl(8)         ! octs 7-8, Ni
   ni = kgds(2)
   kgds(3)=igdstmpl(9)         ! octs 9-10, Nj
   nj = kgds(3)
   kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of 1st grid point
   kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of 1st grid point

   kgds(6)=0                   ! oct 17, resolution and component flags
   if (igdstmpl(1)==2 ) kgds(6)=64
   if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
   if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

   kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.)  ! octs 18-20, Lat of last grid point
   kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.)  ! octs 21-23, Lon of last grid point
   kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.)  ! octs 24-25, di
   kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, dj

   kgds(11) = 0              ! oct 28, scan mode
   if (btest(igdstmpl(19),7)) kgds(11) = 128
   if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
   if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

   kgds(12)=0      ! octs 29-32, reserved
   kgds(19)=0      ! oct 4, # vert coordinate parameters
   kgds(20)=255    ! oct 5, used for thinned grids, set to 255

   res = float(kgds(9)) / 1000.0 * 111.0

 elseif (igdtnum.eq.40) then       !  Gaussian Lat/Lon grid

   iscale=igdstmpl(10)*igdstmpl(11)
   if (iscale==0) iscale=1e6
   kgds(1)=4                   ! oct 6
   kgds(2)=igdstmpl(8)         ! octs 7-8, Ni
   ni = kgds(2)
   kgds(3)=igdstmpl(9)         ! octs 9-10, Nj
   nj = kgds(3)
   kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of 1st grid point
   kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of 1st grid point

   kgds(6)=0                   ! oct 17, resolution and component flags
   if (igdstmpl(1)==2 ) kgds(6)=64
   if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
   if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

   kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point
   kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point
   kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, Di
   kgds(10)=igdstmpl(18)                                 ! octs 26-27, Number of parallels

   kgds(11) = 0              ! oct 28, scan mode
   if (btest(igdstmpl(19),7)) kgds(11) = 128
   if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
   if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

   kgds(12)=0      ! octs 29-32, reserved
   kgds(19)=0      ! oct 4, # vert coordinate parameters
   kgds(20)=255    ! oct 5, used for thinned grids, set to 255

   res = float(kgds(9)) / 1000.0 * 111.0

 elseif (igdtnum.eq.20) then       ! Polar Stereographic Grid

   iscale=1e6
   kgds(1)=5                      ! oct 6, data representation type, polar
   kgds(2)=igdstmpl(8)            ! octs 7-8, nx 
   ni = kgds(2)
   kgds(3)=igdstmpl(9)            ! octs 8-10, ny
   nj = kgds(3)
   kgds(4)=nint(float(igdstmpl(10))/float(iscale)*1000.)  ! octs 11-13, lat of 1st grid point
   kgds(5)=nint(float(igdstmpl(11))/float(iscale)*1000.)  ! octs 14-16, lon of 1st grid point

   kgds(6)=0                      ! oct 17, resolution and component flags
   if (igdstmpl(1) >= 2 .or. igdstmpl(1) <= 5) kgds(6)=64
   if (igdstmpl(1) == 7) kgds(6)=64
   if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128
   if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8

   kgds(7)=nint(float(igdstmpl(14))/float(iscale)*1000.)  ! octs 18-20, lon of orientation
   kgds(8)=nint(float(igdstmpl(15))/float(iscale)*1000.)  ! octs 21-23, dx
   kgds(9)=nint(float(igdstmpl(16))/float(iscale)*1000.)  ! octs 24-26, dy

   kgds(10)=0                ! oct 27, projection center flag
   if (btest(igdstmpl(17),1)) kgds(10) = 128

   kgds(11) = 0              ! oct 28, scan mode
   if (btest(igdstmpl(18),7)) kgds(11) = 128
   if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) +  64
   if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) +  32

   kgds(19)=0    ! oct 4, # vert coordinate parameters
   kgds(20)=255  ! oct 5, used for thinned grids, set to 255

   res = 0.5 * float(kgds(8)+kgds(9)) / 1000. 

 elseif (igdtnum.eq.1) then    ! Rotated Lat/Lon grid

   if (btest(igdstmpl(19),2)) then  ! e-stagger, bit 6 of scan mode is '1'

     iscale=igdstmpl(10)*igdstmpl(11)
     if (iscale == 0) iscale = 1e6
     kgds(1)=203                    ! oct 6, "E" grid
     kgds(2)=igdstmpl(8)            ! octs 7-8, Ni
     ni = kgds(2)
     kgds(3)=igdstmpl(9)            ! octs 9-10, Nj
     nj = kgds(3)
     kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of 1st grid point
     kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of 1st grid point

     kgds(6)=0                      ! oct 17, resolution and component flags
     if (igdstmpl(1)==2 ) kgds(6)=64
     if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
     if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

     kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000  ! octs 18-20, Lat of cent of rotation
     kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.)        ! octs 21-23, Lon of cent of rotation
     kgds(9)=nint(float(igdstmpl(17))/float(iscale)*500.)         ! octs 24-25, Di
                                                                  ! Note!! grib 2 convention twice grib 1
     kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.)       ! octs 26-27, Dj

     kgds(11) = 0                   ! oct 28, scan mode
     if (btest(igdstmpl(19),7)) kgds(11) = 128
     if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
     if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

     kgds(12)=0    ! octs 29-32, reserved
     kgds(19)=0    ! oct 4, # vert coordinate parameters
     kgds(20)=255  ! oct 5, used for thinned grids, set to 255

     res = sqrt( (float(kgds(9)) / 1000.0)**2   +    &
                 (float(kgds(10)) / 1000.0)**2  )
     res = res * 111.0

   else   ! b-stagger

     iscale=igdstmpl(10)*igdstmpl(11)
     if (iscale == 0) iscale = 1e6
     kgds(1)=205                    ! oct 6,     rotated lat/lon for Non-E Stagger grid
     kgds(2)=igdstmpl(8)            ! octs 7-8,  Ni
     ni = kgds(2)
     kgds(3)=igdstmpl(9)            ! octs 9-10, Nj
     nj = kgds(3)
     kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.)  ! octs 11-13, Lat of 1st grid point
     kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.)  ! octs 14-16, Lon of 1st grid point

     kgds(6)=0                      ! oct 17, resolution and component flags
     if (igdstmpl(1)==2 ) kgds(6)=64
     if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) )  kgds(6)=kgds(6)+128
     if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8

     kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.)+90000 ! octs 18-20, Lat of cent of rotation
     kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.)       ! octs 21-23, Lon of cent of rotation
     kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.)       ! octs 24-25, Di
     kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.)      ! octs 26-27, Dj

     kgds(11) = 0                   ! oct 28, scan mode
     if (btest(igdstmpl(19),7)) kgds(11) = 128
     if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) +  64
     if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) +  32

     kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of last grid point
     kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of last grid point

     kgds(19)=0    ! oct 4, # vert coordinate parameters
     kgds(20)=255  ! oct 5, used for thinned grids, set to 255

     res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
             * 0.5 * 111.0

   endif

 else

   print*,'- FATAL ERROR CONVERTING TO GRIB2 GDT.'
   print*,'- UNRECOGNIZED GRID TYPE.'
   call w3tage('SNOW2MDL')
   call errexit(50)

 endif
 
 end subroutine gdt_to_gds

!> Determine length of grib2 gds template array, which is a function of
!! the map projection.
!!      
!! @note call this routine before init_grib2.  
!!
!! @param[in] kgds grib1 gds array
!! @param[in] igdstmplen length of gds template array.
!!
!! condition codes:
!!   47 - unrecognized grid type; fatal
!!
!! @author George Gayno org: w/np2 @date 2014-Sep-28
 subroutine grib2_check (kgds, igdstmplen)
 implicit none

 integer, intent(in)         :: kgds(200)
 integer, intent(  out)      :: igdstmplen

 select case (kgds(1))
   case(4)         ! gaussian
     igdstmplen = 19
   case(203, 205)  ! rotated lat/lon "B" or "E" stagger
     igdstmplen = 22
   case default
     print*,'- FATAL ERROR IN ROUTINE GRIB2_CHECK.'
     print*,'- UNRECOGNIZED GRID TYPE.'
     call w3tage('SNOW2MDL')
     call errexit(47)
 end select

 end subroutine grib2_check

!> Initialize grib2 arrays required by the ncep g2 library according to
!! grib1 gds information. The grib1 gds is held in the kgds array, which
!! is used by the ncep ipolates and w3nco (grib 1) libraries.
!!
!! Call routine grib2_check first to determine igdstmplen.
!!
!! @param[in] century current date/time info
!! @param[in] year current date/time info
!! @param[in] month current date/time info
!! @param[in] day current date/time info
!! @param[in] hour current date/time info
!! @param[in] kgds grib1 gds information
!! @param[in] igdstmplen length of grib2 gdt template.
!! @param[in] lat11 lat of first grid point
!! @param[in] lon11 lon of first grid point
!! @param[in] latlast lat of last grid point
!! @param[in] lonlast lon of last grid point
!! @param[out] igds grib2 section 3 information.
!! @param[out] listsec0 grib2 section 0 information.
!! @param[out] listsec1 grib2 section 1 information.
!! @param[out] ipdsnum grib2 pds template number
!! @param[out] ipdstmpl grib2 pds template array
!! @param[out] igdstmpl grib2 gds template array
!! @param[out] idefnum information for non-reg grid, grid points in each row. 
!! @param[out] ideflist information for non-reg grid, grid points in each row. 
!! @param[out] ngrdpts number of model grid points.
!! @author George Gayno org: w/np2 @date 2014-Sep-28
 subroutine init_grib2(century, year, month, day, hour, kgds, &
                       lat11, latlast, lon11, lonlast, &
                       listsec0, listsec1, igds, ipdstmpl, ipdsnum, igdstmpl,  &
                       igdstmplen, idefnum, ideflist, ngrdpts)
 implicit none

 integer, intent(in   )        :: century, year, month, day, hour
 integer, intent(in   )        :: kgds(200), igdstmplen
 integer, intent(  out)        :: igds(5)
 integer, intent(  out)        :: listsec0(2)
 integer, intent(  out)        :: listsec1(13)
 integer, intent(  out)        :: ipdstmpl(15), ipdsnum
 integer, intent(  out)        :: igdstmpl(igdstmplen)
 integer, intent(  out)        :: idefnum, ideflist
 integer, intent(  out)        :: ngrdpts

 real,    intent(in   )        :: lat11, latlast, lon11, lonlast
 real                          :: scale

! Section 0

 listsec0(1)=0  ! discipline, meteorological fields
 listsec0(2)=2  ! grib version 2

! Section 1

 listsec1(1)=7    ! id of center (ncep)
 listsec1(2)=4    ! subcenter (emc)
 listsec1(3)=8    ! master table version number. wgrib2 does not recognize later tables
 listsec1(4)= 0   ! local table not used
 listsec1(5)= 0   ! signif of ref time - analysis
 if (year == 100) then
   listsec1(6)=century*100 + year
 else
   listsec1(6)=(century-1)*100 + year
 endif
 listsec1(7)=month
 listsec1(8)=day
 listsec1(9)=hour
 listsec1(10:11)=0      ! minutes/secs
 listsec1(12)=0         ! production status of data - ops products
 listsec1(13)=0         ! type of processed products - analysis

! Section 2 - not used

! Section 3 - grid description section

 if (kgds(1) == 4) then  ! gaussian

   igdstmpl(1)=5          ! oct 15; shape of the earth, wgs84
   igdstmpl(2)=255        ! oct 16; scale factor of radius of spherical earth, not used.
   igdstmpl(3)=-1         ! octs 17-20; scale value of radius of spherical earth, not used.
   igdstmpl(4)=255        ! oct 21; scale factor of major axis of elliptical earth, not used.
   igdstmpl(5)=-1         ! octs 22-25; scaled value of major axis of elliptical earth, not used.
   igdstmpl(6)=255        ! oct 26; scale factor of minor axis of elliptical earth, not used.
   igdstmpl(7)=-1         ! octs 27-30; scaled value of minor axis of elliptical earth, not used.
   igdstmpl(8)=kgds(2)    ! octs 31-34; # "i" points
   igdstmpl(9)=kgds(3)    ! octs 35-38; # "j" points
   igdstmpl(10)=1         ! octs 39-42; basic angle
   igdstmpl(11)=10**6     ! octs 43-46; subdivisions of basic angle

   scale=float(igdstmpl(10)*igdstmpl(11))

   igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point

   if (lon11 < 0) then
     igdstmpl(13)=nint((lon11+360.)*scale)  ! octs 51-54; lon of first grid point
   else
     igdstmpl(13)=nint(lon11*scale)
   endif

   igdstmpl(14) = 0  ! oct 55; resolution and component flags
   if (btest(kgds(6),7)) igdstmpl(14) = 48
   if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8

   igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point

   if (lonlast < 0) then
     igdstmpl(16)=nint((lonlast+360.)*scale)  ! octs 60-63; lon of last grid point
   else
     igdstmpl(16)=nint(lonlast*scale)
   endif

   igdstmpl(17)= nint(360.0/float(kgds(2)-1)*scale) ! octs 64-67; di of grid
   igdstmpl(18)= kgds(3)/2 ! octs 68-71; # grid pts between pole and equator

   igdstmpl(19)=0   ! oct 72; scanning mode flag
   if(btest(kgds(11),7)) igdstmpl(19)=128
   if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64
   if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32

   igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1
   igds(2) = kgds(2)*kgds(3)  ! num grid points
   igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid)
   igds(4) = 0 ! regular grid, no appended list
   igds(5) = 40 ! gaussian

   ngrdpts = igds(2)

!  These variables used for non-regular grids.  We are using regular grids
!  (igds(3) equals 0).

   idefnum=1
   ideflist=0

 elseif (kgds(1) == 203 .or. kgds(1) == 205) then

   igdstmpl(1)=5          ! oct 15; shape of the earth, wgs84
   igdstmpl(2)=255        ! oct 16; scale factor of radius of spherical earth, not used.
   igdstmpl(3)=-1         ! octs 17-20; scale value of radius of spherical earth, not used.
   igdstmpl(4)=255        ! oct 21; scale factor of major axis of elliptical earth, not used.
   igdstmpl(5)=-1         ! octs 22-25; scaled value of major axis of elliptical earth, not used.
   igdstmpl(6)=255        ! oct 26; scale factor of minor axis of elliptical earth, not used.
   igdstmpl(7)=-1         ! octs 27-30; scaled value of minor axis of elliptical earth, not used.
   igdstmpl(8)=kgds(2)    ! octs 31-34; # "i" points
   igdstmpl(9)=kgds(3)    ! octs 35-38; # "j" points
   igdstmpl(10)=1         ! octs 39-42; basic angle
   igdstmpl(11)=10**6     ! octs 43-46; subdivisions of basic angle

   scale=float(igdstmpl(10)*igdstmpl(11))

   igdstmpl(12)=nint(lat11*scale) ! octs 47-50; lat of first grid point

   if (lon11 < 0) then
     igdstmpl(13)=nint((lon11+360.)*scale)  ! octs 51-54; lon of first grid point
   else
     igdstmpl(13)=nint(lon11*scale)
   endif

   igdstmpl(14) = 0  ! oct 55; resolution and component flags
   if (btest(kgds(6),7)) igdstmpl(14) = 48
   if (btest(kgds(6),3)) igdstmpl(14) = igdstmpl(14) + 8

   igdstmpl(15)= nint(latlast*scale) ! octs 56-59; lat of last grid point

   if (lonlast < 0) then
     igdstmpl(16)=nint((lonlast+360.)*scale)  ! octs 60-63; lon of last grid point
   else
     igdstmpl(16)=nint(lonlast*scale)
   endif

   if (kgds(1) == 203) igdstmpl(17)= nint(float(kgds(9))*scale/500.)  ! octs 64-67; di of grid.
                                                                      ! iplib "e" grid convention
                                                                      ! is 1/2 the grib convention.
   if (kgds(1) == 205) igdstmpl(17)= nint(float(kgds(9))*scale/1000.) ! octs 64-67; di of grid

   igdstmpl(18)= nint(float(kgds(10))*scale/1000.) ! octs 68-71; dj of grid

   igdstmpl(19)=0   ! oct 72; scanning mode flag
   if(btest(kgds(11),7)) igdstmpl(19)=128
   if(btest(kgds(11),6)) igdstmpl(19)=igdstmpl(19) + 64
   if(btest(kgds(11),5)) igdstmpl(19)=igdstmpl(19) + 32
   if (kgds(1) == 203) igdstmpl(19)=igdstmpl(19) + 4

   igdstmpl(20) = nint(float(kgds(7)-90000)*scale/1000.)  ! octs 73-76; lat of south pole of projection

   if (kgds(8) < 0) then
     igdstmpl(21) = nint(float(kgds(8)+360000)*scale/1000.) ! octs 77-80; long of southern pole of projection.
   else
     igdstmpl(21) = nint(float(kgds(8))*scale/1000.) ! octs 77-80; long of southern pole of projection.
   endif

   igdstmpl(22)=0  ! octs 81-84; angle of rotation of projection

   igds(1) = 0 ! oct 6; source of grid def. specif in table 3.1
   igds(2) = kgds(2)*kgds(3)  ! num grid points
   igds(3) = 0 ! # octets for additional grid pt def (use '0' for regular grid)
   igds(4) = 0 ! regular grid, no appended list
   igds(5) = 1 ! rotated lat/lon

   ngrdpts = igds(2)

!  These variables used for non-regular grids.  We are using regular grids
!  (igds(3) equals 0).

   idefnum=1
   ideflist=0

 end if

! Section 4 - product definition section

 ipdsnum = 0            ! pds template number - table 4.0

 ipdstmpl(1)= 1         ! oct 10; parameter category
! note!! to use a parmeter number >= 192 you must set the local table to '1'
 ipdstmpl(2)= 42        ! oct 11; parameter
 ipdstmpl(3)= 0         ! oct 12; type of generating process
 ipdstmpl(4)= 255       ! oct 13; background generating process identifier
 ipdstmpl(5)= 84        ! oct 14; analysis generating process identifier
 ipdstmpl(6)= 0         ! octs 15-16; hours after ob cutoff
 ipdstmpl(7)= 0         ! oct 17; minutes after ob cutoff
 ipdstmpl(8)= 1         ! oct 18; unit of time range
 ipdstmpl(9)= 0         ! octs 19-22; forecast time in units defined by oct 18
 ipdstmpl(10)=1         ! oct 23; type of first fixed surface
 ipdstmpl(11)=0         ! oct 24; scale factor of first fixed surface
 ipdstmpl(12)=0         ! octs 25-28; scale value of first fixed surface
 ipdstmpl(13)=255       ! oct 29; type of second fixed surface
 ipdstmpl(14)=255       ! oct 30; scale factor of second fixed surface
 ipdstmpl(15)=-2147483647 ! octs 31-34; scaled value of second fixed surface
                              ! note! for these particular octets, using -1 as
                              ! missing does not work because -1 may be an actual
                              ! scaled value.  after looking thru the g2 library
                              ! and some trial and error, i determined that missing
                              ! is minus 2**31-1.

 end subroutine init_grib2

!> Nullify the grib2 gribfield pointers.
!!
!! @param[in] gfld a gribfield data structure
!!
!! @author George Gayno org: w/np2 @date 2014-Sep-28
 subroutine grib2_null(gfld)

 use grib_mod

 implicit none

 type(gribfield), intent(inout)           :: gfld

 nullify(gfld%idsect)
 nullify(gfld%local)
 nullify(gfld%list_opt)
 nullify(gfld%igdtmpl)
 nullify(gfld%ipdtmpl)
 nullify(gfld%coord_list)
 nullify(gfld%idrtmpl)
 nullify(gfld%bmap)
 nullify(gfld%fld)

 end subroutine grib2_null

!> Deallocate the grib2 gribfield pointers.
!!
!! @param[in] gfld a gribfield data structure
!!
!! @author George Gayno org: w/np2 @date 2014-Sep-28
 subroutine grib2_free(gfld)

 use grib_mod

 implicit none

 type(gribfield), intent(inout)    :: gfld

 if (associated(gfld%idsect)) deallocate(gfld%idsect)
 if (associated(gfld%local)) deallocate(gfld%local)
 if (associated(gfld%list_opt)) deallocate(gfld%list_opt)
 if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl)
 if (associated(gfld%ipdtmpl)) deallocate(gfld%ipdtmpl)
 if (associated(gfld%coord_list)) deallocate(gfld%coord_list)
 if (associated(gfld%idrtmpl)) deallocate(gfld%idrtmpl)
 if (associated(gfld%bmap)) deallocate(gfld%bmap)
 if (associated(gfld%fld)) deallocate(gfld%fld)

 end subroutine grib2_free