subroutine grib2_wrt_g2func(fld,bmap,vdate,acc,iptable,lugb,ierr)
!-----------------------------------------------------------------
! ABSTRACT: This routine is to write out a new grib2 file
!   March 2013:     J.Wang
!   23 September 2013: Revised by Youlong Xia for NLDAS grib2 output   
!   21 March 2016: Xia Y.,Revised for NDAS forcing output
!   30 Oct 2017: Y Lin modified for precip analysis 
!     This program is hardwired to write out an 8km cmorph file:
!     grid spec is fixed.  It can be called to write out either the hourly 
!     total, or the sum of 6h total (use the iacc argument)
! Arguments:
!   fld      the (APCP) field
!   bmap     bitmap of fld
!   vdate    valid time yyyymmddhh in character*10
!   acc     accumulating hours (1.0,6.0 etc.)
!   iptable  GRIB2 look-up table to APCP
!   lugb     unit for the grib2 output
!   ierr     error code
!-----------------------------------------------------------------
!
      implicit none
      integer, parameter   :: max_bytes=100000000
!                                       ----:----
      integer, parameter   :: nx=4948,ny=1649
!
      integer iptable(13) ! ipds table
      character vdate*10
      real acc
!
! valid time (idat), reference time (idat0: starting time of accumulation)
!   idat(1)=yyyy
!   idat(2)=mm
!   idat(3)=dd
!   idat(5)=hh
!
      integer idat(8), idat0(8)
! Use w3movdat (rinc(2)=-acc) to compute idat0 from idat. 
      real    rinc(5)
!
      integer listsec0(2)
      integer listsec1(13)
      integer igds(5)
      integer igdstmpllen
      integer ipdstmpllen
      integer idrstmpllen
      integer idrsnum,ibmap,numcoord,ipdsnum,idefnum    
 
      integer,dimension(100) :: igdstmpl
      integer,dimension(100) :: ipdstmpl
      integer,dimension(100) :: idrstmpl
!
      integer ideflist(1)
      real(4) coordlist(1)
!
      character*1 cgrib(max_bytes)
      logical*1 bmap(nx,ny)
!
      real(4),dimension(nx*ny) :: fld
      integer ifilw,i,j,lengrib,lon1,lon2,lat1,lat2,idx,idy,ierr
      integer yyyy,mm,dd,hh,lugb

      real(4) :: dx, dy, alat1, alat2, alon1, alon2
      alon1=0.036378335
      alat1=-59.963614
      dx=0.072756669
      dy=0.072771377
      alon2=alon1+(nx-1)*dx
      alat2=alat1+(ny-1)*dy
!
! code start
!-----------------------------------------------------------------
!
!-- set file unit
      write(6,*) 'grib2_wrt check 1'
      ifilw=lugb
!
      idat=0
      rinc=0.
      read(vdate,'(i4,3i2)') idat(1), idat(2), idat(3), idat(5)
      rinc(2)=-acc
      call w3movdat(rinc,idat,idat0)
      write(6,*) 'idat=', idat
      write(6,*) 'idat0=', idat0

! -------------- add field  -----------------------------------------
      cgrib=''
!
!-- section 0: indicator section 
      listsec0(1)=iptable(1) !Discipline: table 0.0 
!---- (0:Meteorological;1: Hydrlogical; 2:Land) ----------------------
      listsec0(2)=2         ! grib edition number (2:grib2)
!
!-- section 1: identification section
      listsec1(1)=7  ! Identification of orginating center (Table 0)  (7:ncep)
      listsec1(2)=4  ! Identification of orginating subcenter (ON388-Table C) (4:emc)
      listsec1(3)=2         ! GRIB master tables version number (Table 1.0)  (2: Nov 2003)
      listsec1(4)=1         ! Version number of GRIB local tables used to augment Master Tables (Table 1.1)
      listsec1(5)=1         ! Significance of reference time (Table 1.2) (0:ana 1:start of fcst 2:vrfy 3:obs)
      listsec1(6)=idat0(1)  ! Reference time - Year (4 digits)
      listsec1(7)=idat0(2)  ! Reference time - Month
      listsec1(8)=idat0(3)  ! Reference time - Day
      listsec1(9)=idat0(5)  ! Reference time - Hour
      listsec1(10)=0        ! Reference time - Minute
      listsec1(11)=0        ! Reference time - Second
      listsec1(12)=0        ! Production status of data (Table 1.3) (0:opn products 1:opn test products)
      listsec1(13)=1        ! Type of processed data (Table 1.4) (0:ana products 1:fcst products 2:ana & fcst 3: cntl fcst)

      write(6,*) 'grib2_wrt check 3'
       call gribcreate(cgrib,max_bytes,listsec0,listsec1,ierr)
      write(6,*) 'grib2_wrt check 4'
      print*,'gribcreate status=',ierr
!
!-- section 3: grid definition section
      igds(1)=0             ! Source of grid definition (Table 3.0) (0:specified in the code)
      igds(2)=nx*ny         ! Number of grid points in the defined grid
      igds(3)=0             ! Number of octets for optional list of numbers defining number of points 
      igds(4)=0             ! Interpretation of list of numbers defining number of points 
!-- example: Lat/lon
      igds(5)=0            ! Grid definition template number (Table 3.1) (0:Lat/lon)
      igdstmpl=0
      if( igds(5)==0) then
      igdstmpllen=19
!
!-- set up grid definition template 3.0
        igdstmpl=0
        igdstmpl(1)=6       ! Shape of the Earth (Table 3.2) (6:Shape of the Earth = 6,371,229.0 m)
        igdstmpl(8)=nx      ! Ni . number of points along a paralell 
        igdstmpl(9)=ny      ! Nj . number of points along a meridian 
        igdstmpl(10)=0      ! Basic angle of the initial production domain 
        igdstmpl(11)=0      ! Subdivisions of basic angle used to define extreme longitudes and latitudes, and direction increments 
!
        lat1=alat1*1000000
        lon1=alon1*1000000
        lat2=alat2*1000000
        lon2=alon2*1000000
        idx=dx*1000000  
        idy=dy*1000000  

        igdstmpl(12)=lat1   ! La1 - latitude of first grid point
        igdstmpl(13)=lon1   ! Lo1 - longitude of first grid point 
        igdstmpl(14)=48     ! Resolution and component flags (Table 3.3, bits order reversed)
        igdstmpl(15)=lat2   ! La2 - latitude of last grid point
        igdstmpl(16)=lon2   ! Lo2 - longitude of last grid point 
        igdstmpl(17)=idx    ! i direction increment
        igdstmpl(18)=idy    ! j direction increment
        igdstmpl(19)=64     ! Scanning mode (Table 3.4) (+i,+j,i consecutive,row scan same direction)
      endif 
!
      write(6,*) 'grib2_wrt check 5'
      idefnum=1             ! Used if igds(3) .ne. 0. The number of entries in array ideflist
      ideflist=0            ! Used if igds(3) .ne. 0. number of grid points contained in each row ( or column ), Dummy array otherwise
      call addgrid(cgrib,max_bytes,igds,igdstmpl,igdstmpllen,ideflist,
     &             idefnum,ierr)
      write(6,*) 'grib2_wrt check 6'
!      print*,'addgrid status=',ierr
!
!-- section 4: product definition section
      ipdstmpl=0
! ------------ product definition for simultaneous variable ----------
      if(iptable(2).eq.0) then       ! used for simultaneous variables
      ipdsnum=iptable(2)      ! Product Definition Template Number (Table 4.0) 
!(0: at a point in time; 8 for average or accumulation) 
      ipdstmpllen=iptable(3)  ! pdt template length
      ipdstmpl(1)=iptable(4)  ! catogory
      ipdstmpl(2)=iptable(5)  ! parameter
      ipdstmpl(3)=2         ! Type of generating process (Table 4.3) (0:ana, 1:ic, 2:fcst)
      ipdstmpl(4)=0            ! Background generating process identifier 
      ipdstmpl(5)=141          ! Land Data Assimilation and Forecast System identified (ON388TableA) 
      ipdstmpl(6)=0            ! Hours of observational data cutoff after reference time
      ipdstmpl(7)=0            ! Minutes of observational data cutoff after reference time
      ipdstmpl(8)=1            ! Indicator of unit of time range (Table 4.4) (0:minute, 1:hour 2:day)
      ipdstmpl(9)=0            ! Forecast time in units defined by ipdstmpl(8)
      ipdstmpl(10)=iptable(6)  ! Type of first fixed surface (see Code table 4.5) (100:isobaric level)
      ipdstmpl(11)=iptable(7)  ! Scale factor of first fixed surface
      ipdstmpl(12)=iptable(8)  ! Scaled value of first fixed surface
      ipdstmpl(13)=iptable(9)  ! Type of first second surface (see Code table 4.5) (100:isobaric level)
      ipdstmpl(14)=iptable(10) ! Scale factor of second fixed surface
      ipdstmpl(15)=iptable(11) ! Scaled value of second fixed surface
      endif
! ----- product difinition for average or accumulation variable -----
      if(iptable(2).eq.8) then   
      ipdsnum=iptable(2)    ! Product Definition Template
!     Number (Table 4.8) (0: at a point in time; 8 for average or accumulation)
      ipdstmpllen=iptable(3)  ! pdt template length
      ipdstmpl(1)=iptable(4)  ! catogory
      ipdstmpl(2)=iptable(5)  ! parameter
      ipdstmpl(3)=2         ! Type of generating process (Table 4.3)(0:ana, 1:ic, 2:fcst)
      ipdstmpl(4)=0            ! Background generating process identifier
      ipdstmpl(5)=141          ! Land Data Assimilation and Forecast System identified (ON388TableA)
      ipdstmpl(6)=0            ! Hours of observational data cutoff after reference time
      ipdstmpl(7)=0            ! Minutes of observational data cutoff after reference time
      ipdstmpl(8)=1            ! Indicator of unit of time range (Table4.4) (0:minute, 1:hour 2:day)
      ipdstmpl(9)=0            ! Forecast time in units defined by ipdstmpl(8)
      ipdstmpl(10)=iptable(6)  ! Type of first fixed surface (see Code table 4.5) (100:isobaric level)
      ipdstmpl(11)=iptable(7)  ! Scale factor of first fixed surface
      ipdstmpl(12)=iptable(8)  ! Scaled value of first fixed surface
      ipdstmpl(13)=iptable(9)  ! Type of first second surface (see Code table 4.5)
      ipdstmpl(14)=iptable(10) ! Scale factor of second fixed surface
      ipdstmpl(15)=iptable(11) ! Scaled value of second fixed surface
      ipdstmpl(16)=idat(1)     ! End year of overall time interval 
      ipdstmpl(17)=idat(2)     ! End month
      ipdstmpl(18)=idat(3)     ! End day
      ipdstmpl(19)=idat(5)     ! End hour
      ipdstmpl(20)=0           ! End minute
      ipdstmpl(21)=0           ! End second
      ipdstmpl(22)=1           ! number of time ranges
      ipdstmpl(23)=255         ! total number of data values missing
      ipdstmpl(24)=iptable(12) ! average or accumulation
      ipdstmpl(25)=1           ! 1: analysis, 2: forecast from table 4.11
      ipdstmpl(26)=1
      ipdstmpl(27)=1
      ipdstmpl(28)=255
      ipdstmpl(29)=0
      ipdstmpl(30)=0
      endif

      numcoord=0               ! Number of coordinate values after template 
      coordlist=0.             ! Optional list of coordinate values
! ----------- end of Section 4 -----------------------------------------    
!-- section 5: Data Representation Section
      idrstmpl=0
      idrsnum=3             ! Data representation section template number (Table 5.0) (3:Grid Point Data - Complex Packing and Spatial Differencing)
      idrstmpllen=18            ! Length of Data representation section
      idrstmpl(2)=0             ! Binary scale factor
      idrstmpl(3)=iptable(13)   ! Decimal scale factor
      idrstmpl(7)=0             ! Missing value management used (see Code Table 5.5)
      idrstmpl(8)=0             ! Primary missing value substitute
      idrstmpl(9)=0             ! Secondary missing value substitute 
      idrstmpl(17)=1            ! Order of spatial difference (see Code Table 5.6, 2 does not work for Noah output and get segmetation fault) 
!
!-- section 6:       
      ibmap=0  ! Bit-map indicator (Table 6.0) (0:A bit map applies, 255:A bit map doesn't apply)
!
      write(6,*) 'grib2_wrt check 7'
      call addfield(cgrib,max_bytes,ipdsnum,ipdstmpl,ipdstmpllen,
     &              coordlist,numcoord,idrsnum,idrstmpl,idrstmpllen,
     &              fld,nx*ny,ibmap,bmap,ierr)
      if(ierr.ne.0) then
      print*,'addfield fails and status=',ierr
      endif
!      print*,'addfield status=',ierr

!-- finalize  GRIB message after all section
!-- adds the End Section ( "7777" )
      call gribend(cgrib,max_bytes,lengrib,ierr)
!      print*,'gribend status=',ierr
!      print*,'length of the final GRIB2 message in octets =',lengrib
!
      call wryte(ifilw, lengrib, cgrib)

      end subroutine grib2_wrt_g2func