!
!  Uses getgb2 to unpack records
!
      use grib_mod
      type(gribfield) :: gfldper,gfldwht
      parameter (maxpts=1000000,lcgrib=4*maxpts)
      integer,dimension(200) :: jids,jpdt,jgdt
c.. Fields for encodeing into grib2
      CHARACTER(len=1) :: cgrib(lcgrib),cin(lcgrib)
      integer :: listsec0(3)
      integer :: listsec1(13)
      integer :: igds(5),igdstmpl(200),ipdstmpl(200)
      integer :: ideflist,idefnum
      integer :: idrstmpl(200)
      real :: coordlist(1)
      real :: wstpn(maxpts)
      Logical*1 :: bmap(maxpts)
      logical :: unpack=.true.
      character*11 envvar
      character(len=80) :: gfile,gfilei,g2file
      real vtest,pi,gr,const1
!
c
      lugb = 10   !  gribfile
      lugi = 11   !  grib index file
!!!
      envvar='XLFUNIT_   '
      write(envvar(9:10),fmt='(I2)') lugb
      call getenv(envvar,gfile)
      call baopenr(lugb,gfile,iret1)
c
      envvar='XLFUNIT_   '
      write(envvar(9:10),fmt='(I2)') lugi
      call getenv(envvar,gfilei)
      call baopenr(lugi,gfilei,iret2)
c
!     Output Grib2 file
c
      IFL2=20   ! output grib2 gridded file
      envvar='XLFUNIT_   '
      write(envvar(9:10),fmt='(I2)') IFL2
      call getenv(envvar,g2file)
      CALL BAOPENW(ifl2,g2file,IOS)
c
!    Extract Wave height & Period ONLY!!
c
       do ifcts = 1,1
c
c----------------------------------------------------
!  Get Wave Height
c
         jskp = 0
         jdisc=-1
         jpdtn=0
         jgdtn=-1
         jids=-9999
         jpdt=-9999
         jgdt=-9999
c
        jpdt(1) = 0
        jpdt(2) = 3
        jpdt(10) = 1

      call getgb2(lugb,lugi,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,
     &            unpack,jskp,gfldwht,iret)
         if ( iret.ne.0) then
            print *,' getgb2 whts = ',iret
            cycle
         endif
c
           fldmax=-9999999.
           fldmin=9999999.
           sum=0.0
            idd = 0
c
              do j=1,gfldwht%ngrdpts
                if (.not. gfldwht%bmap(j)) cycle
                if (gfldwht%fld(j).gt.fldmax) fldmax=gfldwht%fld(j)
                if (gfldwht%fld(j).lt.fldmin) fldmin=gfldwht%fld(j)
                sum=sum+gfldwht%fld(j)
                 idd = idd + 1
              enddo
c
      print*,'Data Values SWhgt:',' dd',idd,' ndpts=',gfldwht%ndpts
           write(6,fmt='("MIN=",f21.8," AVE=",f21.8,
     &          " MAX=",f21.8)') fldmin,sum/gfldwht%ndpts,fldmax
c
!  Get Wave Period      
c
         jskp = 0
         jdisc=-1
         jpdtn=0
         jgdtn=-1
         jids=-9999
         jpdt=-9999
         jgdt=-9999
c
        jpdt(1) = 0
        jpdt(2) = 11
        jpdt(10) = 1

      call getgb2(lugb,lugi,jskp,jdisc,jids,jpdtn,jpdt,jgdtn,jgdt,
     &            unpack,jskp,gfldper,iret)
         if ( iret.ne.0) then
            print *,' getgb2 wper = ',iret
            cycle
         endif
!
           fldmax=-9999999.
           fldmin=9999999.
           sum=0.0
             idd = 0
c
              do j=1,gfldper%ngrdpts
                if (.not. gfldper%bmap(j)) cycle
                if (gfldper%fld(j).gt.fldmax) fldmax=gfldper%fld(j)
                if (gfldper%fld(j).lt.fldmin) fldmin=gfldper%fld(j)
                sum=sum+gfldper%fld(j)
                idd = idd + 1
              enddo
c
       print*,'Data Values Wper:',' dd',idd,' ndpts=',gfldper%ndpts,
     &        'ngrdpts= ',gfldper%ngrdpts
           write(6,fmt='("MIN=",f21.8," AVE=",f21.8,
     &          " MAX=",f21.8)') fldmin,sum/gfldper%ndpts,fldmax
c
c----------------------------------------------------
c        Compute Wave Steepness
c set up some constant values to compute wave steepness
c
c
            pi=3.141592654
            gr=9.802
            const1=gr/(2*pi)
c
         do  jj = 1,gfldper%ngrdpts
          if (gfldper%fld(jj).gt. 0.0 .and.
     *       gfldper%fld(jj).lt. 100.)then
             const2=1./(const1*gfldper%fld(jj)*gfldper%fld(jj)) 
             vtest = gfldwht%fld(jj) * const2
              wstpn(jj) = vtest
              bmap(jj) = gfldper%bmap(jj)
c
            else
               wstpn(jj) = gfldper%fld(jj)
                bmap(jj) = gfldper%bmap(jj)  
c
          endif
         enddo
c
       print*,' Array Wstpness Completed'
c
!   Replace parameter number with wsteep
c
             gfldper%ipdtmpl(2) = 192
             
c
! Pack record into new grib2 file
c
c
!  create grib message
c
!==    Initialize new GRIB2 message and pack
! GRIB2 sections 0 (Indicator Section) and 1 (Identification  Section)
!==
c
        listsec0 = 0
c
       listsec0(1)= 10 ! Discipline-GRIB Master Table Number (see Code Table 0.0)
c                       ! Met = 0 & Ocean = 10
       listsec0(2)=2 ! GRIB Edition Number (currently 2)
       listsec0(3)=0
c
       listsec1(1)=7       ! 7  Id of orginating centre (Common Code Table C-1)
       listsec1(2)=0 !"EMC"! Id of orginating sub-centre (local table)/Table C/on388
       listsec1(3)=2       ! GRIB Master Tables Version Number (Code Table 1.0)
       listsec1(4)=1       !per Brent! GRIB Local Tables Version Number(Code Table 1.)
       listsec1(5)=1       ! Significance of Reference Time (Code Table 1.2)
       listsec1(6)= gfldper%idsect(6)    ! Reference Time - Year (4 digits)
       listsec1(7)= gfldper%idsect(7)    ! Reference Time - Month
       listsec1(8)= gfldper%idsect(8)    ! Reference Time - Day
       listsec1(9)= gfldper%idsect(9)    ! Reference Time - Hour(cycle:0,6,12,18)
       listsec1(10)= 0      ! Reference Time - Minute
       listsec1(11)= 0      ! Reference Time - Second

       listsec1(12)=0      ! Production status of data (Code Table 1.3)
       listsec1(13)=1      ! Type of processed data (Code Table 1.4)
c
      call gribcreate(cgrib,lcgrib,listsec0,listsec1,ierr)
      if (ierr.ne.0) then
        write(6,*) ' ERROR creating new GRIB2 field = ',ierr
        stop 2
      endif
c
! assume predefined grid for simplicity
!
!==> Pack up Grid Definition Section (Section 3) add to GRIB2 message.
c
c    Lat/long grid
c
                igds = 0
c
           igds(1)=0                  !Source of grid definition (see Code Table 3.0)
           igds(2)=gfldper%ngrdpts    ! num of grid points
           igds(3)=0   !Number of octets needed for each additional grid points dfn
           igds(4)=00  !Interp of list for optional points defn (Code Table 3.11)
           igds(5)=0          !Grid Definition Template Number (Code Table 3.1)
c
           idefnum = 0
           ideflist=0       !Used if igds(3) .ne. 0. Dummy array otherwise
c
            igdstmpl = 00
c
           igdstmpl(1)=6
           igdstmpl(2)=0
           igdstmpl(3)=0
           igdstmpl(4)=0
           igdstmpl(5)=0
           igdstmpl(6)=0
           igdstmpl(7)=0
           igdstmpl(8)=gfldper%igdtmpl(8)    ! num points along parallel
           igdstmpl(9)=gfldper%igdtmpl(9)    ! num points along meridian
           igdstmpl(10)=0
           igdstmpl(11)=0
           igdstmpl(12)=gfldper%igdtmpl(12)   ! lat of first grid point
           igdstmpl(13)=gfldper%igdtmpl(13)   ! convert W to E
cccc
           igdstmpl(14)=gfldper%igdtmpl(14)    ! Resolution and Component flags
           igdstmpl(15)=gfldper%igdtmpl(15)  ! lagdtmpl(j)% of last grid point
ccccc
           igdstmpl(16)=gfldper%igdtmpl(16)    ! convert W to E
cc
           igdstmpl(17)=gfldper%igdtmpl(17)   ! Increment of lat
           igdstmpl(18)=gfldper%igdtmpl(18)  ! Increment of long
           igdstmpl(19)=gfldper%igdtmpl(19)       ! scanning mode
cccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      call addgrid(cgrib,lcgrib,igds,igdstmpl,200,ideflist,
     &             idefnum,ierr)
      if (ierr.ne.0) then
        write(6,*) ' ERROR adding GRIB2 grid = ',ierr
        stop 3
      endif
c
! assume PDT 4.0 for simplicity
c
       ipdsnum=0    !  anly or fcst at a hozontal level or level
c
! Parameter category (see Code Table 4.1)
c
         ipdstmpl = 00
c
        ipdstmpl(1) = 0   ! 0 = ocean / 2 = meteo
c
c    Get parameter number (see Code Table 4.2)
c
        ipdstmpl(2) = 192   ! wave steepness
c
c  Type of generating process: analysis or forecast(see code Table 4.3)
c
c anal(0), init(1), fcst(2), bias corr fcst(3), ensemb fcst(4), etc
c
c  For  Model/Forecast fields
c
        ipdstmpl(3)=2     ! forecast
        ipdstmpl(4)=0     !background generating process identifier
!                         (defined by originating Center)
c
        ipdstmpl(5)= 11 ! Generating process or Model Number
!                            (defined by originating Center)
c
        ipdstmpl(6)=0     ! hours of observational data cutoff after reference time
        ipdstmpl(7)=0     ! minutes of observational data cutoff after reference time
c
        ipdstmpl(8)=gfldper%ipdtmpl(8)  !indic of unit of time range( Code Table 4.4)
        ipdstmpl(9)= gfldper%ipdtmpl(9)  !forecast time in units defined by pdstmpl(8)
c
c
        ipdstmpl(10)=gfldper%ipdtmpl(10)  ! type of level(see Code Table 4.5)1st level
        ipdstmpl(11)=gfldper%ipdtmpl(11)        ! scale factor of pdstmpl(10)
        ipdstmpl(12)=gfldper%ipdtmpl(12)        ! scaled value of pdstmpl(10)
c
        ipdstmpl(13)=255   ! type of level (See Code Table 4.5) 2nd level
        ipdstmpl(14)=0     ! scale factor of ipdstmpl(13)
        ipdstmpl(15)=0     ! scaled value of ipdstmpl(13)
!
ccccccccccccccccccccc
      numcoord=0
        coordlist= 0.0     !needed for hybrid vertical coordinate
! set bitmap flag
          ibmap=gfldper%ibmap    ! Bitmap indicator ( see Code Table 6.0 )
      idrsnum=40    ! Data Rep Template Number ( see Code Table 5.40 ; grd pnt data)
c
      idrstmpl=0
c
!********************************************************************************
! idrstmpl(1): reference value(R) (IEEE 32-bit floating-point value)     
! idrstmpl(2): binary scale factor (E)                                           *
! idrstmpl(3): decimal scale factor (D)                                          *
! idrstmpl(4): number of bits used for each packed value for simple packing      *
!              or for each group reference value for complex packing or          *
!              spatial differencing                                              *
! idrstmpl(5): type of original field values (See Code Table 5.1)                *
!********************************************************************************
c
      idrstmpl(3)=gfldper%idrtmpl(3) ! binary scale factor (E)
      idrstmpl(4)=gfldper%idrtmpl(4)
c
         ndpts = gfldper%ngrdpts
c
      call addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,200,
     &              coordlist,numcoord,idrsnum,idrstmpl,200,
     &              wstpn,ndpts,ibmap,bmap,ierr)
      if (ierr.ne.0) then
          write(6,*) ' ERROR adding GRIB2 field = ',ierr
          stop 4
      endif
c
c
       print*,'IPDSTMPL ',ipdsnum,idrsnum,ipdstmpl(1:25)
c
       print*,'IGDSTMPL ',igds(1:5),igdstmpl(1:22)
c
!
! End GRIB2 field
!
      call gribend(cgrib,lcgrib,lengrib,ierr)
      if (ierr.ne.0) then
        write(6,*) ' ERROR ending new GRIB2 message = ',ierr
        stop 5
      endif
      print *,' writing ',lengrib,' bytes...'
       call wryte(ifl2,lengrib,cgrib)
c

! Free memory when done with field
c
        call gf_free(gfldper)
        call gf_free(gfldwht)
c
C  Get next forecast for wave hgt and period 
c
      Enddo       !   get next fcst
c
        CALL BACLOSE (LUGB,iret)
         CALL BACLOSE (LUGI,jret)
          CALL BACLOSE(ifl2,kret)
      print*,'iret jret kret ',iret,jret,kret
c
      stop
      end