!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE GRIBDUMP                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE gribdump(nx,ny,nz,nstyps, nchanl,filnam, grdbas,             & 1,106
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,              &
           ubar,vbar,ptbar,pbar,rhobar,qvbar,                           &
           x,y,z,zp,hterain, j1,j2,j3,                                  &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           tem1,item1,temxy1,temxy2,itemxy1,itemxy2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write history data into channel nchanl as GRIB data.
!
!  All data read in are located at the original staggered grid points.
!
!  Note: coordinate fields are dumped as 3 dimensional fields which
!  have been converted from meters to kilometers.  This is for the
!  convenience of the plotting applications.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/01/1995
!
!  MODIFICATION HISTORY:
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  1/07/2000 (Eric Kemp)
!  Snow cover is now written as GRIB variable 66 (snow depth (m))
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    nchanl   FORTRAN I/O channel number for history data output.
!    filnam   File name to store the GRIB messages
!    grdbas   Flag indicating if this is a call for the data dump
!             of grid and base state arrays only. If so, grdbas=1.
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Vertical component of Cartesian velocity at a given
!             time level (m/s)
!    ptprt    Perturbation potential temperature at a given time
!             level (K)
!    pprt     Perturbation pressure at a given time level (Pascal)
!    qv       Water vapor specific humidity at a given time level (kg/kg)
!    qc       Cloud water mixing ratio at a given time level (kg/kg)
!    qr       Rainwater mixing ratio at a given time level (kg/kg)
!    qi       Cloud ice mixing ratio at a given time level (kg/kg)
!    qs       Snow mixing ratio at a given time level (kg/kg)
!    qh       Hail mixing ratio at a given time level (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    hterain  Terrain height (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Skin temperature at the ground or ocean surface (K)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT:
!
!    None.
!
!  WORK ARRAY:
!
!    tem1     Temporary work array.
!    item1    Integer temporary work array.
!    temxy1   2-D temporary work array.
!    temxy2   2-D temporary work array.
!    itemxy1  2-D integer temporary work array.
!    itemxy2  2-D integer temporary work array.
!
!
!-----------------------------------------------------------------------
!
!  The following parameters are passed into this subroutine through
!  a common block in globcst.inc, and they determine which
!  variables are output.
!
!  grdout =0 or 1. If grdout =0, grid variables are not dumped.
!  basout =0 or 1. If basout =0, base state variables are not dumped.
!  varout =0 or 1. If varout =0, model variables are not dumped.
!  mstout =0 or 1. If mstout =0, water variables are not dumped.
!  rainout=0 or 1. If rainout=0, rain variables are not dumped.
!  prcout =0 or 1. If prcout =0, precipitation rates are not dumped.
!  iceout =0 or 1. If iceout =0, qi, qs and qh are not dumped.
!  tkeout =0 or 1. If tkeout =0, tke is not dumped.
!  trbout =0 or 1. If trbout =0, kmh and kmv are not dumped.
!  sfcout =0 or 1. If sfcout =0, surface variables are not dumped.
!  landout=0 or 1. If landout=0, surface propertty arrays are not dumped.
!  radout =0 or 1. If radout =0, precipitation rates are not dumped.
!  flxout =0 or 1. If flxout =0, precipitation rates are not dumped.
!
!  These following parameters are also passed in through common
!  blocks in globcst.inc.
!
!  runname,curtim,umove,vmove,xgrdorg,ygrdorg
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INCLUDE 'mp.inc'          ! Message passing parameters.

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  INTEGER :: nchanl            ! FORTRAN I/O channel number for output
  CHARACTER (LEN=*     ) :: filnam ! File name
  INTEGER :: grdbas            ! If this is a grid/base state array dump

  REAL :: u     (nx,ny,nz)     ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz)     ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz)     ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)

  REAL :: qv    (nx,ny,nz)     ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz)     ! Turbulent Kinetic Energy ((m/s)**2)

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  REAL :: hterain(nx,ny)       ! Terrain height.

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as - d( zp )/d( x )
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as - d( zp )/d( y )
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z )

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)  !Soil type fraction
  INTEGER :: vegtyp (nx,ny)          ! Vegetation type
  REAL :: lai    (nx,ny)          ! Leaf Area Index
  REAL :: roufns (nx,ny)          ! Surface roughness
  REAL :: veg    (nx,ny)          ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)    ! Temperature at surface (K)
                                     ! (in top 1 cm layer)
  REAL :: tsoil  (nx,ny,0:nstyps)    ! Deep soil temperature (K)
                                     ! (in deep 1 m layer)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)    ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount
  REAL :: snowdpth(nx,ny)            ! Snow depth (m)

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rates (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1(nx,ny,nz)       ! Temporary work array
  REAL :: temxy1(nx,ny)        ! temporary work array
  REAL :: temxy2(nx,ny)        ! temporary work array

  INTEGER :: item1(nx,ny,nz)   ! Integer temporary work array
  INTEGER :: itemxy1(nx,ny)    ! Integer temporary work array
  INTEGER :: itemxy2(nx,ny)    ! Integer temporary work array
!
!-----------------------------------------------------------------------
!
!  Parameters for GRIB packing
!
!-----------------------------------------------------------------------
!
  INTEGER :: nbufsz
  PARAMETER   ( nbufsz = 800000 )  ! Size of GRIB buffer

  INTEGER :: ipds(25)          ! PDS integer array
  INTEGER :: igdss(25)         ! GDS integer array for s-grid
  INTEGER :: igdsu(25)         ! GDS integer array for u-grid
  INTEGER :: igdsv(25)         ! GDS integer array for v-grid

  INTEGER :: ibdshd(4)         ! BDS header

  CHARACTER (LEN=1) :: pds(28)       ! PDS ( GRIB Section 1)
  CHARACTER (LEN=1) :: gds(42)       ! GDS ( GRIB Section 2)
  CHARACTER (LEN=1) :: bds(nbufsz)   ! BDS ( GRIB Section 4)
  INTEGER :: ibds(nbufsz/4)    ! identical to BDS

  EQUIVALENCE ( bds,ibds )

  CHARACTER (LEN=1) :: mgrib(nbufsz) ! Buffer to carry GRIB messages

  INTEGER :: nbytes            ! Byte position to start read or write
  INTEGER :: wdlen             ! Length of machine word
!
!-----------------------------------------------------------------------
!
!  Parameters describing this routine
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=40) :: fmtver
  PARAMETER (fmtver='004.10 ARPS Data Format 10')
!
!-----------------------------------------------------------------------
!
!  Variables to determine the length of machine word, wdlen
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=8) :: ctema,ctemb

  INTEGER :: itema,itemb

  EQUIVALENCE ( itema,ctema )
  EQUIVALENCE ( itemb,ctemb )
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,is
  INTEGER :: itype     ! Array type: 0 - floating, 1 - integer
  INTEGER :: dflag     ! Dimension flag: 0 - 3-D, 1 - 2-D
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!  include 'gribcst.inc'
!
!-----------------------------------------------------------------------
!
!  Save ipds & igds for future use
!
!-----------------------------------------------------------------------
!
  LOGICAL :: firstcall(0:mgrdmax)

  SAVE firstcall, wdlen
  DATA firstcall/ mgrdmax*.true.,.true. /
  DATA ctema / '12345678' /
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!   Initialize the integer array for BDS header
!
!-----------------------------------------------------------------------
!
  ibdshd(1) = 0
  ibdshd(2) = 0
  ibdshd(3) = 0
  ibdshd(4) = 0
!
!-----------------------------------------------------------------------
!
!  Initialize the integer array, ipds, for PDS.
!
!-----------------------------------------------------------------------
!
  CALL mkipds( nx,ny,nz, ipds )
!
!-----------------------------------------------------------------------
!
!  Initialize the GDS integer arrays for different grid type, igdsu,
!  igdsv, and igdss
!
!-----------------------------------------------------------------------
!
  CALL mkigds( nx,ny,nz, 0, igdss )
  CALL mkigds( nx,ny,nz, 1, igdsu )
  CALL mkigds( nx,ny,nz, 2, igdsv )

  IF (mgrid < 0 .or. mgrid > mgrdmax) THEN
    mgrid = 0
  ENDIF

  IF ( firstcall(mgrid) ) THEN
!
!-----------------------------------------------------------------------
!
!  Calculate the length of machine word, wdlen, in bytes, which will
!  be used to pack the data. Assume the length has only two possible
!  value: 4 and 8
!
!-----------------------------------------------------------------------
!
    itemb = itema
    IF ( ctema /= ctemb ) THEN
      wdlen = 4
    ELSE
      wdlen = 8
    END IF
    IF (myproc == 0) WRITE (6,'(a,i2,a)')  &
        'The length of machine word is ',wdlen,' bytes'
!
!-----------------------------------------------------------------------
!
!  Check if number of bits per datum fits size of computer word.
!
!-----------------------------------------------------------------------
!
    IF ( wdlen*8 < grbpkbit ) THEN
      IF (myproc == 0) WRITE (6,'(a,i3,a,i3/a,a/a)')  &
          'Number of bits per datum, ',grbpkbit,                        &
          ', exceeds the machine word length, ',wdlen*8,                &
          'Please change grbpkbit in input namelist &output ', &
          'to fit the machine word length.',                            &
          'Program stoped in GRIBDUMP'
      CALL arpsstop(" ",1)
    END IF
!
!-----------------------------------------------------------------------
!
!  Write the sample of GrADS control file used to display the ARPS
!  GRIB files. This control file is to display the variables on the
!  scalar points. The GRIB grid type was given by ipds(5)
!
!-----------------------------------------------------------------------
!
    CALL gribcntl( nx,ny,nz,x,y,z,zp, ipds(5), temxy1,temxy2 )

    firstcall(mgrid) = .false.

  END IF

  IF (myproc == 0) WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim
!
!-----------------------------------------------------------------------
!
!  Write the header of ARPS GRIB file
!
!-----------------------------------------------------------------------
!
  CALL wrthishd(nx,ny,nz,nchanl,grdbas,fmtver, x,y,z,                   &
                wdlen,nbufsz,mgrib,nbytes)
!
!-----------------------------------------------------------------------
!
!  If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
  IF(grdout == 1 .OR. grdbas == 1 ) THEN

    itype = 0       ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D array
    ipds(8)  = 8    ! geometric height (m)
    ipds(9)  = 103  ! z-coorinates
    ipds(18) = 0    ! forecast time = 0 for grid and base variables

    IF (myproc == 0) WRITE(6,'(a)') 'Writing zp'
    CALL edgfill(zp,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
    CALL wrtgrib(nx,ny,nz,1,z,zp,item1,                                 &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF    ! grdout
!
!-----------------------------------------------------------------------
!
!  If basout=1, write out base state variables.
!
!-----------------------------------------------------------------------
!
  IF(basout == 1 .OR. grdbas == 1 ) THEN

    itype = 0       ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D array
    ipds(9)  = 103  ! z-coorinates
    ipds(18) = 0    ! forecast time = 0 for grid and base variables

    ipds(8)  = 33   ! u wind (m/s)
    IF (myproc == 0) WRITE(6,'(a,f10.5)') 'Writing ubar ',ubar(1,1,1)
    CALL edgfill(ubar,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,ubar,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdsu,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 34   ! v wind (m/s)
    IF (myproc == 0) WRITE(6,'(a,f10.5)') 'Writing vbar ',vbar(1,1,1)
    CALL edgfill(vbar,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,vbar,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdsv,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)


    ipds(8)  = 13   ! potential temperature (K)
    IF (myproc == 0) WRITE(6,'(a,f10.2)') 'Writing ptbar ',ptbar(1,1,1)
    CALL edgfill(ptbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,ptbar,item1,                              &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 1    ! pressure (Pascal)
    IF (myproc == 0) WRITE(6,'(a,f10.2)') 'Writing pbar ', pbar(1,1,1)
    CALL edgfill(pbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,pbar,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    IF(mstout == 1) THEN

      ipds(8)  = 51   ! specific humidity (kg/kg)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing qvbar'
      CALL edgfill(qvbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL wrtgrib(nx,ny,nz,0,z,qvbar,item1,                            &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

    END IF

    IF(landout == 1) THEN

      dflag = 1         ! 2-D array
      ipds(9)  = 111    ! depth below land surface

      IF( nstyp <= 1 ) THEN

        itype = 1       ! integer array
        ibdshd(3) = itype

        ipds(11) = 0    ! 0 centimeter for surface

        ipds(8)  = 224  ! soil type
        IF (myproc == 0) WRITE(6,'(a)') 'Writing soiltyp '

        CALL iedgfill(soiltyp(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,tem1(1,1,1),soiltyp(1,1,1),            &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ELSE

        itype = 1       ! integer array
        ibdshd(3) = itype
        ipds(8)  = 224  ! soil type

        IF (myproc == 0) WRITE(6,'(a)') 'Writing soiltyp '

        DO is=1,nstyp
          ipds(11) = is ! 0 centimeter for surface

          CALL iedgfill(soiltyp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,       &
                        1,1,1,1)
          CALL wrtgrib(nx,ny,1,0,z,tem1(1,1,1),soiltyp(1,1,is),         &
                       nchanl,nbytes,wdlen,itype,dflag,                 &
                       ipds,igdss,ibdshd,pds,gds,                       &
                       nbufsz,bds,ibds,                                 &
                       mgrib,temxy1,temxy2,itemxy1,itemxy2)
        END DO

        itype = 0       ! floating array
        ibdshd(3) = itype
        ipds(8)  = 226  ! redefine 226 for fraction of soil types

        IF (myproc == 0) WRITE(6,'(a)') 'Writing stypfrct'

        DO is=1,nstyp
          ipds(11) = is ! 0 centimeter for surface

          CALL edgfill(stypfrct(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,       &
                       1,1,1,1)
          CALL wrtgrib(nx,ny,1,0,z,stypfrct(1,1,is),item1(1,1,1),       &
                       nchanl,nbytes,wdlen,itype,dflag,                 &
                       ipds,igdss,ibdshd,pds,gds,                       &
                       nbufsz,bds,ibds,                                 &
                       mgrib,temxy1,temxy2,itemxy1,itemxy2)
        END DO

      END IF

      itype = 1       ! integer array
      ibdshd(3) = itype

      ipds(11) = 0
      ipds(8)  = 225  ! vegtyp
      IF (myproc == 0) WRITE(6,'(a)') 'Writing vegtyp'
      CALL iedgfill(vegtyp ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,tem1(1,1,1),vegtyp,                      &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      itype = 0       ! floating array
      ibdshd(3) = itype

      ipds(8)  = 227  ! redefine 227 for LAI
      IF (myproc == 0) WRITE(6,'(a)') 'Writing lai'
      CALL edgfill(lai    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,lai,item1(1,1,1),                        &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 83   ! surface roughness (m)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing roufns'
      CALL edgfill(roufns ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,roufns,item1(1,1,1),                     &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 87   ! vegetation (%)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing veg'
      CALL edgfill(veg    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      DO j = 1,ny
        DO i = 1,nx
          tem1(i,j,1) = veg(i,j) * 100.
        END DO
      END DO
      CALL wrtgrib(nx,ny,1,0,z,tem1(1,1,1),item1(1,1,1),                &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

    END IF

  END IF

  IF ( grdbas == 1 ) RETURN

  ipds(18) = nint(curtim/60) ! forecast time in minutes
!
!-----------------------------------------------------------------------
!
!  If varout = 1, Write out u, v, w, pt, and p.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Write out u,v and w.
!
!-----------------------------------------------------------------------
!
  IF(varout == 1) THEN

    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0                  ! 3-D array
    ipds(9)  = 103             ! z-coorinates

    ipds(8)  = 33   ! u wind (m/s)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing u'
    CALL edgfill(u,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,u,item1,                                  &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdsu,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 34   ! v wind (m/s)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing v'
    CALL edgfill(v,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,v,item1,                                  &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdsv,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 40   ! w wind (m/s)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing w'
    CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
    CALL wrtgrib(nx,ny,nz,1,z,w,item1,                                  &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
    ipds(8)  = 13   ! potential temperature (K)
    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tem1(i,j,k) = ptbar(i,j,k) + ptprt(i,j,k)
        END DO
      END DO
    END DO
    IF (myproc == 0) WRITE(6,'(a,f10.2)') 'Writing pt ', tem1(1,1,1)
    CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,tem1,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 1    ! pressure (Pascal)
    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tem1(i,j,k) = pbar(i,j,k) + pprt(i,j,k)
        END DO
      END DO
    END DO
    CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    IF (myproc == 0) WRITE(6,'(a,f10.2)') 'Writing p ', tem1(1,1,1)
    CALL wrtgrib(nx,ny,nz,0,z,tem1,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF     ! varout
!
!-----------------------------------------------------------------------
!
!  If mstout = 1, write out moisture scalars.
!
!-----------------------------------------------------------------------
!
  IF(mstout == 1) THEN

    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0                  ! 3-D array
    ipds(9)  = 103             ! z-coorinates

    ipds(8)  = 51   ! specific humidity (kg/kg)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing qv'
    CALL edgfill(qv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,qv,item1,                                 &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 153  ! cloud water mixing ratio (kg/kg)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing qc'
    CALL edgfill(qc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,qc,item1,                                 &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 170  ! rain water mixing ratio (kg/kg)
    IF (myproc == 0) WRITE(6,'(a)') 'Writing qr'
    CALL edgfill(qr,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,qr,item1,                                 &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    IF(iceout == 1) THEN

      dflag = 0      ! 3-D
      ipds(9)  = 103  ! z-coorinates

      ipds(8)  = 178  ! ice mixing ratio (kg/kg)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing qi'
      CALL edgfill(qi,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL wrtgrib(nx,ny,nz,0,z,qi,item1,                               &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 171  ! snow mixing ratio (kg/kg)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing qs'
      CALL edgfill(qs,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL wrtgrib(nx,ny,nz,0,z,qs,item1,                               &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 179  ! graupel hail mixing ratio (kg/kg)
      IF (myproc == 0) WRITE(6,'(a)') 'Writing qh'
      CALL edgfill(qh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL wrtgrib(nx,ny,nz,0,z,qh,item1,                               &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

    END IF   !iceout

    IF(rainout == 1) THEN

      dflag = 1       ! 2-D
      ipds(9)  = 1    ! at surface

      ipds(8)  = 62   ! large scale precipitation for raing (kg/m**2)
      ipds(11) = 0    ! at surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing raing'
      CALL edgfill(raing,   1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,raing,item1,                             &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 63   ! convective precipitation for rainc (kg/m**2)
      ipds(11) = 0    ! at surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing rainc'
      CALL edgfill(rainc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,rainc,item1,                             &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

    END IF   !rainout

    IF ( prcout == 1 ) THEN

      CALL edgfill(prcrate, 1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)

      dflag = 1       ! 2-D
      ipds(8)  = 59   ! precipitation rates (kg/m**2/s)
      ipds(9)  = 111  ! use 111 to save more than one precip. rate

      ipds(11) = 1    ! total precipitation rate
      IF (myproc == 0) WRITE(6,'(a)') 'Writing prcrate1'
      CALL wrtgrib(nx,ny,1,0,z,prcrate(1,1,1),item1,                    &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(11)  = 2   ! grid scale precipitation rates
      IF (myproc == 0) WRITE(6,'(a)') 'Writing prcrate2'
      CALL wrtgrib(nx,ny,1,0,z,prcrate(1,1,2),item1,                    &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(11) = 3   ! cumulus precipitation rates
      IF (myproc == 0) WRITE(6,'(a)') 'Writing prcrate3'
      CALL wrtgrib(nx,ny,1,0,z,prcrate(1,1,3),item1,                    &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(11) = 4   ! microphysics precipitation rates
      IF (myproc == 0) WRITE(6,'(a)') 'Writing prcrate4'
      CALL wrtgrib(nx,ny,1,0,z,prcrate(1,1,4),item1,                    &
                    nchanl,nbytes,wdlen,itype,dflag,                    &
                    ipds,igdss,ibdshd,pds,gds,                          &
                    nbufsz,bds,ibds,                                    &
                    mgrib,temxy1,temxy2,itemxy1,itemxy2)

    END IF   ! prcout

  END IF   !mstout

  IF( tkeout == 1 ) THEN
    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates
    ipds(8)  = 158  ! tke (J kg**-1)
    CALL edgfill(tke,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,tke,item1,                                &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF
!
!-----------------------------------------------------------------------
!
!  If trbout = 1, write out the turbulence parameter, km.
!
!-----------------------------------------------------------------------
!
  IF( trbout == 1 ) THEN
    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates
    ipds(8)  = 185  ! redefine 185 for kmh
    IF (myproc == 0) WRITE(6,'(a)') 'Writing kmh'
    CALL edgfill(kmh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,kmh,item1,                                &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates
    ipds(8)  = 186  ! redefine 186 for kmv
    CALL edgfill(kmv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,kmv,item1,                                &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF   ! trbout

!
!-----------------------------------------------------------------------
!
!  If sfcout = 1, write out the surface variables,
!  tsfc, tsoil, wetsfc, wetdp, and wetcanp.
!
!-----------------------------------------------------------------------
!
  IF ( sfcout == 1 ) THEN
    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 1         ! 2-D
    ipds(9)  = 111    ! below land surface

    IF( nstyp <= 1 ) THEN

      ipds(8)  = 85   ! soil temperature (K)
      ipds(11) = 1    ! 1 cm below the surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing tsfc'
      CALL edgfill(tsfc(1,1,0),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,tsfc(1,1,0),item1,                       &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 85   ! soil temperature (K)
      ipds(11) = 100  ! 100 cm below the surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing tsoil'
      CALL edgfill(tsoil(1,1,0),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,tsoil(1,1,0),item1,                      &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
      ipds(11) = 1    ! 1 cm below the surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing wetsfc'
      CALL edgfill(wetsfc(1,1,0),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,wetsfc(1,1,0),item1,                     &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
      ipds(11) = 100  ! 100 cm below the surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing wetdp'
      CALL edgfill(wetdp(1,1,0),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,wetdp(1,1,0),item1,                      &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)

      ipds(8)  = 223  ! amount of liquid water retained on canopy (m)
      ipds(11) = 0    ! 0 cm below the surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing wetcnpy'
      CALL edgfill(wetcanp(1,1,0),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,wetcanp(1,1,0),item1,                    &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ELSE

      DO is=0,nstyp
        ipds(8)  = 85   ! soil temperature (K)
        ipds(11) = 1+is ! 1 cm below the surface
        IF (myproc == 0) WRITE(6,'(a)') 'Writing tsfc'
        CALL edgfill(tsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,              &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,tsfc(1,1,is),item1,                    &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)

        ipds(8)  = 85   ! soil temperature (K)
        ipds(11) = 100+is ! 100 cm below the surface
        IF (myproc == 0) WRITE(6,'(a)') 'Writing tsoil'
        CALL edgfill(tsoil(1,1,is),  1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,tsoil(1,1,is),item1,                   &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)

        ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
        ipds(11) = 1+is ! 1 cm below the surface
        IF (myproc == 0) WRITE(6,'(a)') 'Writing wetsfc'
        CALL edgfill(wetsfc(1,1,is), 1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,wetsfc(1,1,is),item1,                  &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)

        ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
        ipds(11) = 100+is ! 100 cm below the surface
        IF (myproc == 0) WRITE(6,'(a)') 'Writing wetdp'
        CALL edgfill(wetdp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,            &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,wetdp(1,1,is),item1,                   &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)

        ipds(8)  = 223  ! amount of liquid water retained on canopy (m)
        ipds(11) = is   ! 0 cm below the surface
        IF (myproc == 0) WRITE(6,'(a)') 'Writing wetcnpy'
        CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,           &
                     1,1,1,1)
        CALL wrtgrib(nx,ny,1,0,z,wetcanp(1,1,is),item1,                 &
                     nchanl,nbytes,wdlen,itype,dflag,                   &
                     ipds,igdss,ibdshd,pds,gds,                         &
                     nbufsz,bds,ibds,                                   &
                     mgrib,temxy1,temxy2,itemxy1,itemxy2)
      END DO
    END IF

    IF (snowout == 1) THEN
      itype = 0                  ! floating array
      ibdshd(3) = itype
      ipds(8)  = 66   ! Snow depth (m)
      ipds(9)  = 1    ! at surface
      ipds(11) = 0    ! at surface
      IF (myproc == 0) WRITE(6,'(a)') 'Writing snowdpth'
      CALL edgfill(snowdpth,1,nx,1,nx-1,1,ny,1,ny-1, 1,1,1,1)
      CALL wrtgrib(nx,ny,1,0,z,snowdpth,item1,                          &
                   nchanl,nbytes,wdlen,itype,dflag,                     &
                   ipds,igdss,ibdshd,pds,gds,                           &
                   nbufsz,bds,ibds,                                     &
                   mgrib,temxy1,temxy2,itemxy1,itemxy2)
    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  If radout = 1, write out radiation arrays
!
!-----------------------------------------------------------------------
!
  IF( radout == 1 ) THEN
    itype = 0       ! integer array
    ibdshd(3) = itype

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates
    ipds(8)  = 216  ! radiation forcing, radfrc
    CALL edgfill(radfrc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL wrtgrib(nx,ny,nz,0,z,radfrc,item1,                             &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    dflag = 1       ! 2-D
    ipds(9)  = 1    ! at surface
    ipds(8)  = 204  ! downward short wave radiation flux
    ipds(11) = 0    ! 0 cm below the surface
    CALL edgfill(radsw,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,radsw,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    dflag = 1       ! 2-D
    ipds(9)  = 1    ! at surface
    ipds(8)  = 232  ! downward total (net) radiation flux
    CALL edgfill(rnflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,rnflx,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF
!
!-----------------------------------------------------------------------
!
!  If flxout = 1, write out surface fluxes
!
!-----------------------------------------------------------------------
!
  IF( flxout == 1 ) THEN
    itype = 0       ! integer array
    ibdshd(3) = itype

    dflag = 1       ! 2-D
    ipds(9)  = 1    ! at surface
    ipds(8)  = 124  ! u flux
    ipds(11) = 0    ! 0 cm below the surface
    CALL edgfill(usflx,1,nx,1,nx, 1,ny,1,ny-1, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,usflx,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 125  ! v flux
    CALL edgfill(vsflx,1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,vsflx,item1,                               &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 122  ! sensible heat flux for ptsflx
    CALL edgfill(ptsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,ptsflx,item1,                              &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

    ipds(8)  = 121  ! latent heat flux for qvsflx
    CALL edgfill(qvsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL wrtgrib(nx,ny,1,0,z,qvsflx,item1,                              &
                  nchanl,nbytes,wdlen,itype,dflag,                      &
                  ipds,igdss,ibdshd,pds,gds,                            &
                  nbufsz,bds,ibds,                                      &
                  mgrib,temxy1,temxy2,itemxy1,itemxy2)

  END IF

  RETURN
END SUBROUTINE gribdump
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GRIBREAD                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE gribread(nx,ny,nz,nstyps, grdbas, inch, time, x,y,z,zp,      & 2,63
           uprt, vprt, wprt, ptprt, pprt,                               &
           qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                      &
           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,                &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in GRIB data set created by ARPS using history dump format
!  No. 10.
!  All data read in are located at the original staggered grid points
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/01/1995
!
!  MODIFICATION HISTORY:
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  1/07/2000 (Eric Kemp)
!  Subroutine now reads in GRIB variable 66 (snow depth (m))
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    grdbas   Data read flag.
!             =1, only grid and base state arrays will be read
!             =0, all arrays will be read based on data
!                 parameter setting.
!    inch     Channel number for binary reading.
!             This channel must be opened for unformatted reading
!             by the calling routine.
!
!  OUTPUT:
!
!    time     Time in seconds of data read from "filename"
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    uprt     x component of perturbation velocity (m/s)
!    vprt     y component of perturbation velocity (m/s)
!    wprt     Vertical component of perturbation velocity in
!             Cartesian coordinates (m/s).
!    ptprt    Perturbation potential temperature (K)
!    pprt     Perturbation pressure (Pascal)
!
!    qvprt    Perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state air density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    ireturn  Return status indicator
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!  TEMPORARY:
!
!    tem1     Temporary work array.
!    item1    Integer temporary work array.
!    temxy1   2-D temporary work array.
!    itemxy1  2-D integer temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: x(nx)                ! x-coord. of the physical and compu
                               ! -tational grid. Defined at u-point(m).
  REAL :: y(ny)                ! y-coord. of the physical and compu
                               ! -tational grid. Defined at v-point(m).
  REAL :: z(nz)                ! z-coord. of the computational grid.
                               ! Defined at w-point on the staggered
                               ! grid(m).
  REAL :: zp(nx,ny,nz)         ! Physical height coordinate defined at
                               ! w-point of the staggered grid(m).

  INTEGER :: grdbas            ! Data read flag.
  INTEGER :: inch              ! Channel number for binary reading
  REAL :: time                 ! Time in seconds of data read
                               ! from "filename"

  REAL :: uprt  (nx,ny,nz)     ! Perturbation u-velocity (m/s)
  REAL :: vprt  (nx,ny,nz)     ! Perturbation v-velocity (m/s)
  REAL :: wprt  (nx,ny,nz)     ! Perturbation w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)

  REAL :: qvprt (nx,ny,nz)     ! Perturbation water vapor mixing
                               ! ratio (kg/kg)
  REAL :: qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
  REAL :: tke   (nx,ny,nz)     ! Turbulent Kinetic Energy ((m/s)**2)
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: wbar  (nx,ny,nz)     ! Base state w-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor mixing ratio

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)  ! Soil type fraction
  INTEGER :: vegtyp (nx,ny)          ! Vegetation type
  REAL :: lai    (nx,ny)          ! Leaf Area Index
  REAL :: roufns (nx,ny)          ! Surface roughness
  REAL :: veg    (nx,ny)          ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)    ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyps)    ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)    ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount
  REAL :: snowdpth(nx,ny)            ! Snow depth (m)

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rates (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL, allocatable :: tem1(:,:,:)     ! Temporary work array
  REAL, allocatable :: temxy1(:,:)     ! temporary work array

  INTEGER, allocatable :: item1(:,:,:) ! Integer temporary work array
  INTEGER, allocatable :: itemxy1(:,:) ! Integer temporary work array

  INTEGER :: ireturn
!
!-----------------------------------------------------------------------
!
!  Parameters describing routine that wrote the gridded data
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Parameters for GRIB packing
!
!-----------------------------------------------------------------------
!
  INTEGER :: nbufsz
  PARAMETER ( nbufsz = 800000 )  ! Size of GRIB buffer

  INTEGER :: ipds(25)          ! PDS integer array
  INTEGER :: igdss(25)         ! GDS integer array for s-grid
  INTEGER :: igdsu(25)         ! GDS integer array for u-grid
  INTEGER :: igdsv(25)         ! GDS integer array for v-grid
  INTEGER :: ibdshd(4)         ! Integer array for BDS header

  CHARACTER (LEN=1) :: bds(nbufsz)   ! BDS ( GRIB Section 4)
  INTEGER :: ibds(nbufsz/4)    ! identical S ( GRIB Section 4)

  EQUIVALENCE ( bds,ibds )

  CHARACTER (LEN=1) :: mgrib(nbufsz) ! Buffer to carry GRIB messages

  INTEGER :: itype             ! Data type indicator
  INTEGER :: dflag             ! Dimension flag
!
!-----------------------------------------------------------------------
!
!  Variables to determine the length of machine word, wdlen
!
!-----------------------------------------------------------------------
!
  INTEGER :: wdlen             ! Length of machine word

  CHARACTER (LEN=8) :: ctema,ctemb

  INTEGER :: itema,itemb

  EQUIVALENCE ( itema,ctema )
  EQUIVALENCE ( itemb,ctemb )

  LOGICAL :: firstcall

  SAVE firstcall, wdlen
  DATA firstcall/ .true. /
  DATA ctema / '12345678' /
!
!-----------------------------------------------------------------------
!
!  Parameters describing this routine
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=40) :: fmtver
  PARAMETER (fmtver='004.10 ARPS Data Format 10')
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: lchanl
  PARAMETER (lchanl=6)      ! Channel number for formatted printing.

  INTEGER :: i,j,k,is

  REAL :: tema1, tema2, temb1, temb2
  REAL :: alatnot(2)

  INTEGER :: istat, nstyp1, nstyp2
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE (tem1(nx,ny,nz),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "GRIBREAD: ERROR allocating tem1, returning"
    RETURN
  END IF
  ALLOCATE (item1(nx,ny,nz),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "GRIBREAD: ERROR allocating item1, returning"
    RETURN
  END IF
  ALLOCATE (temxy1(nx,ny),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "GRIBREAD: ERROR allocating temxy1, returning"
    RETURN
  END IF
  ALLOCATE (itemxy1(nx,ny),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "GRIBREAD: ERROR allocating itemxy1, returning"
    RETURN
  END IF

  IF ( firstcall ) THEN

    itemb = itema
    IF ( ctema /= ctemb ) THEN
      wdlen = 4
    ELSE
      wdlen = 8
    END IF

    WRITE (6,'(a,i2,a)')                                                &
        'The machine word length is ',wdlen,' bytes'

    firstcall = .false.

  END IF
!
!-----------------------------------------------------------------------
!
!  Check header info
!
!-----------------------------------------------------------------------
!
  nstyp2 = nstyp

  CALL rdhishd( nx,ny,nz,inch,grdbas,fmtver,time,x,y,z,                 &
                nbufsz,mgrib )

  nstyp1 = nstyp ! value read in from data
  nstyp = nstyp2 ! original value
  IF ( nstyp1 < 1 ) THEN
    nstyp1 = 1
  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate the length of machine word, wdlen, in bytes, which will
!  be used to pack the data. Assume the length has only two possible
!  value: 4 and 8
!
!-----------------------------------------------------------------------
!
  dx = x(2) - x(1)
  dy = y(2) - y(1)

  alatnot(1) = trulat1
  alatnot(2) = trulat2

  CALL setmapr(mapproj, 1.0, alatnot, trulon)
  CALL lltoxy( 1,1, ctrlat,ctrlon, tema1, tema2 )

  temb1 = tema1 - (FLOAT(nx-3)/2.) * dx
  temb2 = tema2 - (FLOAT(ny-3)/2.) * dy

  CALL setorig( 1, temb1, temb2)

  CALL setcornerll( nx,ny, x,y )               ! set corner lat/lon
!
!-----------------------------------------------------------------------
!
!   Initialize the integer array for BDS header
!
!-----------------------------------------------------------------------
!
  ibdshd(1) = 0
  ibdshd(2) = 0
  ibdshd(3) = 0
  ibdshd(4) = 0
!
!-----------------------------------------------------------------------
!
!  Get the integer array, ipds, from the header.
!
!-----------------------------------------------------------------------
!
  CALL mkipds( nx,ny,nz, ipds )
!
!-----------------------------------------------------------------------
!
!  Get the GDS integer arrays for different grid type, igdsu, igdsv,
!  and igdss, from the header
!
!-----------------------------------------------------------------------
!
  CALL mkigds( nx,ny,nz, 0, igdss )
  CALL mkigds( nx,ny,nz, 1, igdsu )
  CALL mkigds( nx,ny,nz, 2, igdsv )
!
!-----------------------------------------------------------------------
!
!  Get zp
!
!----------------------------------------------------------------------
!
  IF( grdin == 1 .OR. grdbas == 1 ) THEN

    itype = 0       ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D array
    ipds(8)  = 8    ! geometric height (m)
    ipds(9)  = 103  ! z-coorinates
    ipds(18) = 0    ! forecast time = 0 for grid and base variables

    CALL rdgrib(nx,ny,nz,1,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                zp,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' zp.'
  END IF  ! grdin
!
!-----------------------------------------------------------------------
!
!  Read in base state fields
!
!----------------------------------------------------------------------
!
  IF( basin == 1 .OR. grdbas == 1 ) THEN

    itype = 0       ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D array
    ipds(9)  = 103  ! z-coorinates
    ipds(18) = 0    ! forecast time = 0 for grid and base variables

    ipds(8)  = 33   ! u wind (m/s)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdsu,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                ubar,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' ubar.'

    ipds(8)  = 34   ! v wind (m/s)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdsv,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                vbar,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' vbar.'

    ipds(8)  = 13   ! potential temperature (K)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                ptbar,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' ptbar.'

    ipds(8)  = 1    ! pressure (Pascal)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                pbar,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' pbar.'

    IF( mstin == 1) THEN
      ipds(8)  = 51   ! specific humidity (kg/kg)
      CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                       &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  qvbar,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' qvbar.'
    END IF

    IF(landin == 1) THEN

      dflag = 1       ! 2-D array
      ipds(9)  = 111  ! depth below land surface

      IF( nstyp1 <= 1 ) THEN
        itype = 1       ! integer array
        ibdshd(3) = itype

        ipds(8)  = 224  ! soil type
        ipds(11) = 0    ! 0 centimeter for surface

        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    tem1,soiltyp(1,1,1), temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' soiltyp.'

      ELSE

        itype = 1       ! integer array
        ibdshd(3) = itype
        ipds(8)  = 224  ! soil type

        DO is=1,nstyp1

          ipds(11) = is ! number to identify soil types

          CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                    &
                      ipds,igdss,ibdshd,                                &
                      nbufsz,bds,ibds,mgrib,                            &
                      tem1,soiltyp(1,1,is), temxy1,itemxy1)
          WRITE(lchanl,910) ipds(8),' soiltyp.'

        END DO

        itype = 0       ! floating array
        ibdshd(3) = itype
        ipds(8)  = 226  ! redefine 226 for fraction of soil types

        DO is=1,nstyp1

          ipds(11) = is ! fraction of different soil types

          CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                    &
                      ipds,igdss,ibdshd,                                &
                      nbufsz,bds,ibds,mgrib,                            &
                      stypfrct(1,1,is),item1, temxy1,itemxy1)
          WRITE(lchanl,910) ipds(8),' stypfrct.'

        END DO

      END IF

      CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct)

      itype = 1       ! integer array
      ibdshd(3) = itype

      ipds(8)  = 225  ! vegtyp
      ipds(11) = 0
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  tem1,vegtyp, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' vegtyp.'

      itype = 0       ! floating array
      ibdshd(3) = itype

      ipds(8)  = 227  ! redefine 227 for LAI
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  lai,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' lai.'

      ipds(8)  = 83   ! surface roughness (m)
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  roufns,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' roufns.'

      ipds(8)  = 87   ! vegetation (%)
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  tem1,item1, temxy1,itemxy1)

      DO j = 1,ny
        DO i = 1,nx
          veg(i,j) = tem1(i,j,1) / 100.
        END DO
      END DO
      WRITE(lchanl,910) ipds(8),' veg.'

    END IF

  END IF

  IF( grdbas == 1 ) GO TO 930

  ipds(18) = nint(time/60) ! forecast time in minutes

  IF( varin == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!
!----------------------------------------------------------------------
!
    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0                  ! 3-D array
    ipds(9)  = 103             ! z-coorinates

    ipds(8)  = 33   ! u wind (m/s)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdsu,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tem1,item1, temxy1,itemxy1)

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          uprt(i,j,k) = tem1(i,j,k) - ubar(i,j,k)
        END DO
      END DO
    END DO
    WRITE(lchanl,910) ipds(8),' uprt.'

    ipds(8)  = 34   ! v wind (m/s)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdsv,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tem1,item1, temxy1,itemxy1)

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          vprt(i,j,k) = tem1(i,j,k) - vbar(i,j,k)
        END DO
      END DO
    END DO
    WRITE(lchanl,910) ipds(8),' vprt.'

    ipds(8)  = 40   ! w wind (m/s)
    CALL rdgrib(nx,ny,nz,1,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                wprt,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' wprt.'
!
!-----------------------------------------------------------------------
!
!  Read in scalars
!
!----------------------------------------------------------------------
!
    ipds(8)  = 13   ! potential temperature (K)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tem1,item1, temxy1,itemxy1)

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          ptprt(i,j,k) = tem1(i,j,k) - ptbar(i,j,k)
        END DO
      END DO
    END DO
    WRITE(lchanl,910) ipds(8),' ptprt.'

    ipds(8)  = 1    ! pressure (Pascal)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tem1,item1, temxy1,itemxy1)

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          pprt(i,j,k) = tem1(i,j,k) - pbar(i,j,k)
        END DO
      END DO
    END DO
    WRITE(lchanl,910) ipds(8),' pprt.'

  END IF
!
!-----------------------------------------------------------------------
!
!  Read in moisture variables
!
!----------------------------------------------------------------------
!
  IF( mstin == 1 ) THEN

    itype = 0                  ! floating array
    ibdshd(3) = itype

    dflag = 0                  ! 3-D array
    ipds(9)  = 103             ! z-coorinates

    ipds(8)  = 51   ! specific humidity (kg/kg)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tem1,item1, temxy1,itemxy1)

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          qvprt(i,j,k) = tem1(i,j,k) - qvbar(i,j,k)
        END DO
      END DO
    END DO
    WRITE(lchanl,910) ipds(8),' qvprt.'

    ipds(8)  = 153  ! cloud water mixing ratio (kg/kg)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                qc,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' qc.'

    ipds(8)  = 170  ! rain water mixing ratio (kg/kg)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                qr,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' qr.'

    IF( icein == 1 ) THEN

      dflag = 0       ! 3-D
      ipds(9)  = 103  ! z-coorinates

      ipds(8)  = 178  ! ice mixing ratio (kg/kg)
      CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                       &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  qi,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' qi.'

      ipds(8)  = 171  ! snow mixing ratio (kg/kg)
      CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                       &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  qs,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' qs.'

      ipds(8)  = 179  ! graupel hail mixing ratio (kg/kg)
      CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                       &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  qh,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' qh.'

    END IF

    IF( rainin == 1 ) THEN
      dflag = 1       ! 2-D
      ipds(9)  = 1    ! at surface

      ipds(8)  = 62   ! large scale precipitation for raing (kg/m**2)
      ipds(11) = 0    ! at surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  raing,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' raing.'

      ipds(8)  = 63   ! convective precipitation for rainc (kg/m**2)
      ipds(11) = 0    ! at surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  rainc,item1, temxy1,itemxy1)

      WRITE(lchanl,910) ipds(8),' rainc.'
    END IF

    IF( prcin == 1 ) THEN
      dflag = 1       ! 2-D
      ipds(8)  = 59   ! precipitation rates (kg/m**2/s)
      ipds(9)  = 111  ! use 111 for saving more than one precip.

      ipds(11) = 1    ! total precipitation rate
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  prcrate(1,1,1),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' prcrate1.'

      ipds(11) = 2    ! grid scale precipitation rate
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  prcrate(1,1,2),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' prcrate2.'

      ipds(11) = 3    ! cumulus precipitation rate
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  prcrate(1,1,3),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' prcrate3.'

      ipds(11) = 4    ! microphysics precipitation rate
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  prcrate(1,1,4),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' prcrate4.'

    END IF

  END IF

  IF( tkein == 1 ) THEN

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates

    ipds(8)  = 158  ! tke (J kg**-1)
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                tke,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' tke.'

  END IF

  IF( trbin == 1 ) THEN

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates

    ipds(8)  = 185  ! redefine 185 for kmh
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                kmh,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' kmh.'

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates

    ipds(8)  = 186  ! redefine 186 for kmv
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                kmv,item1, temxy1,itemxy1)

    WRITE(lchanl,910) ipds(8),' kmv.'

  END IF

  IF( sfcin == 1 ) THEN

    dflag = 1         ! 2-D
    ipds(9)  = 111    ! below land surface

    IF( nstyp1 <= 1 ) THEN

      ipds(8)  = 85   ! soil temperature (K)
      ipds(11) = 1    ! 1 cm below the surface

      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  tsfc(1,1,0),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' tsfc.'

      ipds(8)  = 85   ! soil temperature (K)
      ipds(11) = 100  ! 100 cm below the surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  tsoil(1,1,0),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' tsoil.'

      ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
      ipds(11) = 1    ! 1 cm below the surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  wetsfc(1,1,0),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' wetsfc.'

      ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
      ipds(11) = 100  ! 100 cm below the surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  wetdp(1,1,0),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' wetdp.'

      ipds(8)  = 223  ! amount of liquid water retained on canopy (m)
      ipds(11) = 0    ! 0 cm below the surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  wetcanp(1,1,0),item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' wetcanp.'

    ELSE

      DO is=0,nstyp1

        ipds(8)  = 85   ! soil temperature (K)
        ipds(11) = 1+is ! 1 cm below the surface
        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    tsfc(1,1,is),item1, temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' tsfc.'

        ipds(8)  = 85   ! soil temperature (K)
        ipds(11) = 100+is ! 100 cm below the surface
        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    tsoil(1,1,is),item1, temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' tsoil.'

        ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
        ipds(11) = 1+is ! 1 cm below the surface
        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    wetsfc(1,1,is),item1, temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' wetsfc.'

        ipds(8)  = 144  ! Soil moisture in (m**3/m**3)
        ipds(11) = 100+is ! 100 cm below the surface
        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    wetdp(1,1,is),item1, temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' wetdp.'

        ipds(8)  = 223  ! amount of liquid water retained on canopy (m)
        ipds(11) = is   ! 0 cm below the surface
        CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                      &
                    ipds,igdss,ibdshd,                                  &
                    nbufsz,bds,ibds,mgrib,                              &
                    wetcanp(1,1,is),item1, temxy1,itemxy1)
        WRITE(lchanl,910) ipds(8),' wetcanp.'

      END DO

    END IF

    CALL fix_soil_nstyp(nx,ny,nstyp1,nstyp,tsfc,tsoil,wetsfc,wetdp,wetcanp)

    IF (snowcin == 1) THEN  ! Snow cover is depricated
      itype = 1
      ibdshd(3) = itype
      ipds(8)  = 143  ! Snow cover
      ipds(9)  = 111
      ipds(11) = 0
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  tem1,item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' snowcvr -- discarding.'
    END IF

    IF (snowin == 1) THEN
      itype = 0       ! floating array
      ibdshd(3) = itype
      ipds(8)  = 66   ! Snow depth (m)
      ipds(9)  = 1    ! at surface
      ipds(11) = 0    ! at surface
      CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                        &
                  ipds,igdss,ibdshd,                                    &
                  nbufsz,bds,ibds,mgrib,                                &
                  snowdpth,item1, temxy1,itemxy1)
      WRITE(lchanl,910) ipds(8),' snowdpth.'
    END IF

  END IF

  IF( radin == 1 ) THEN
    itype = 0       ! floating array
    ibdshd(3) = itype

    dflag = 0       ! 3-D
    ipds(9)  = 103  ! z-coorinates

    ipds(8)  = 216  ! radiation forcing, radfrc
    CALL rdgrib(nx,ny,nz,0,z,dflag,inch, wdlen,                         &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                radfrc,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' radfrc.'

    dflag = 1       ! 2-D
    ipds(9)  = 1    ! at surface
    ipds(11) = 0    ! at surface

    ipds(8)  = 204  ! downward short wave radiation flux
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                radsw,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' radsw.'

    ipds(8)  = 232  ! downward total (net) radiation flux
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                rnflx,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' rnflx.'

  END IF

  IF( flxin == 1 ) THEN
    itype = 0       ! integer array
    ibdshd(3) = itype

    ipds(11) = 0    ! at surface
    ipds(8)  = 124  ! u flux
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                usflx,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' usflx.'

    ipds(8)  = 125  ! v flux
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                vsflx,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' vsflx.'

    ipds(8)  = 122  ! sensible heat flux for ptsflx
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                ptsflx,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' ptsflx.'

    ipds(8)  = 121  ! latent heat flux for qvsflx
    CALL rdgrib(nx,ny,1,0,z,dflag,inch, wdlen,                          &
                ipds,igdss,ibdshd,                                      &
                nbufsz,bds,ibds,mgrib,                                  &
                qvsflx,item1, temxy1,itemxy1)
    WRITE(lchanl,910) ipds(8),' qvsflx.'

  END IF

  910   FORMAT(1X,'Field ',i4,' was read into array',a)

!
!-----------------------------------------------------------------------
!
!  Friendly exit message
!
!----------------------------------------------------------------------
!
  930   CONTINUE

  WRITE(6,'(/a,F8.1,a/)')                                               &
      ' Data at time=', time/60,' (min) were successfully read.'

  ireturn = 0

  DEALLOCATE (tem1,stat=istat)
  DEALLOCATE (item1,stat=istat)
  DEALLOCATE (temxy1,stat=istat)
  DEALLOCATE (itemxy1,stat=istat)

  RETURN

END SUBROUTINE gribread
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE WRTGRIB                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE wrtgrib(nx,ny,nz,wstag,z,fvar,ivar,                          & 51,3
           nchanl,nbytes,wdlen,itype,dflag,                             &
           ipds,igds,ibdshd,pds,gds,                                    &
           nbufsz,bds,ibds,                                             &
           mgrib,tem1,tem2,item1,item2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write variable, fvar for floating array or ivar for integer array,
!  into file pointer nchanl which points to a GRIB file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/01/1995
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       number of grid points in the vertical
!
!    wstag    Staggering flag of w-grid
!    z        z coordinate of grid points in computational space (m)
!
!    fvar     Floating array to be written into GRIB file
!    ivar     Integer  array to be written into GRIB file
!
!    nchanl   FILE pointer of GRIB stream file which was opened by
!             a C program, GOPEN
!
!    itype    Data type: 0 for floating array and 1 for integer array
!    dflag    Dimension flag: 0 for 3-D array and 1 for 2-D array
!    nbytes   Byte position at which the GRIB message starts
!    wdlen    Length of machine word (i.e. size of integer) in bytes
!
!    ipds     Integer array to carry the parameters for generating
!             PDS (Section 1).
!    igds     Integer array to carry the parameters for generating
!             GDS (Section 3).
!
!    pds      PDS (GRIB Section 1)
!    gds      GDS (GRIB Section 3)
!
!    nbufsz   Size of buffer of a GRIB message array
!    bds      BDS (GRIB Section 4)
!    mgrib    Working buffer of GRIB message array
!
!  OUTPUT:
!
!    nbtyes   Position at which the next GRIB message starts
!
!  WORK ARRAY:
!
!    tem1     Temporary work array.
!    item1    Integer temporary work array.
!    tem2     Temporary work array.
!    item2    Integer temporary work array.
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in x,y,z dir.

  INTEGER :: wstag             ! Staggering flag for w-grid

  REAL :: z(nz)                ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: fvar(nx,ny,nz)    ! Floating array
  INTEGER :: ivar(nx,ny,nz)    ! Integer  array

  INTEGER :: nchanl            ! FILE pointer indicates the GRIB file
  INTEGER :: itype             ! Data type: 0 - floating, 1 - integer
  INTEGER :: dflag             ! Dimension flag: 0 - 3-D, 1 - 2-D
  INTEGER :: nbytes            ! Starting byte position
  INTEGER :: wdlen             ! Length of machine word

  INTEGER :: ipds(25)          ! ipds
  INTEGER :: igds(25)          ! igds
  INTEGER :: ibdshd(4)         ! BDS header

  CHARACTER (LEN=1) :: pds(28)       ! PDS
  CHARACTER (LEN=1) :: gds(42)       ! GDS

  INTEGER :: nbufsz            ! Size of GRIB message array
  INTEGER :: ibds(nbufsz/wdlen) ! Identical to BDS
  CHARACTER (LEN=1) :: bds(nbufsz)   ! BDS
  CHARACTER (LEN=1) :: mgrib(nbufsz) ! Working buffer

  REAL :: tem1 (nx,ny)         ! Temporary work array
  REAL :: tem2 (nx,ny)         ! Temporary work array

  INTEGER :: item1(nx,ny)      ! Integer temporary work array
  INTEGER :: item2(nx,ny)      ! Integer temporary work array
!
!-----------------------------------------------------------------------
!
!  GRIB parameters
!
!-----------------------------------------------------------------------
!
  INTEGER :: msglen    ! Length of GRIB message
  INTEGER :: npts
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, wstag1
  REAL :: tema
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  npts = nx*ny

  IF ( dflag == 1 ) THEN

    DO j = 1,ny
      DO i = 1,nx
        tem1 (i,j) = fvar(i,j,1)
        item1(i,j) = ivar(i,j,1)
      END DO
    END DO

    CALL gribenc(itype,wdlen,grbpkbit,npts,tem1,item1,                  &
                 ipds,igds,ibdshd,pds,gds,                              &
                 nbufsz,bds,ibds,msglen,mgrib,                          &
                 tem2,item2)

    WRITE (nchanl) (mgrib(i),i=1,msglen)

    RETURN
  END IF

  IF ( wstag > 1 .OR. wstag < 0 ) THEN
    WRITE(6,'(a,i1)')                                                   &
        'Do not know the w-grid type, wstag = ',wstag,                  &
        'Reset it to scalar grid.'
    wstag1 = 0
  ELSE
    wstag1 = wstag
  END IF

  DO k = 1,nz-1
    DO j = 1,ny
      DO i = 1,nx
        tem1 (i,j) = fvar(i,j,k)
        item1(i,j) = ivar(i,j,k)
      END DO
    END DO

    IF ( wstag1 == 0 ) THEN
      tema = 0.5 * (z(k) + z(k+1))
    ELSE
      tema = z(k)
    END IF
    ipds(11) = nint(tema)

    CALL gribenc(itype,wdlen,grbpkbit,npts,tem1,item1,                  &
                 ipds,igds,ibdshd,pds,gds,                              &
                 nbufsz,bds,ibds,msglen,mgrib,                          &
                 tem2,item2)

    WRITE (nchanl) (mgrib(i),i=1,msglen)

  END DO

  DO j = 1,ny
    DO i = 1,nx
      tem1 (i,j) = fvar(i,j,nz)
      item1(i,j) = ivar(i,j,nz)
    END DO
  END DO

  ipds(11) = nint(
  CALL gribenc(itype,wdlen,grbpkbit,npts,tem1,item1,                    &
               ipds,igds,ibdshd,pds,gds,                                &
               nbufsz,bds,ibds,msglen,mgrib,                            &
               tem2,item2)

  WRITE (nchanl) (mgrib(i),i=1,msglen)

  RETURN
END SUBROUTINE wrtgrib
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RDGRIB                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE rdgrib(nx,ny,nz,wstag,z,dflag,nchanl, wdlen,                 & 52,3
           ipds,igds,ibdshd,                                            &
           nbufsz,bds,ibds,mgrib,                                       &
           fvar,ivar, tem1,item1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write variable, fvar for floating array or ivar for integer array,
!  into file pointer nchanl which points to a GRIB file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/01/1995
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       number of grid points in the vertical
!
!    dflag    Dimension flag, 0: 3-D, 1: 2-D
!    wstag    Staggering flag of w-grid
!
!    nchanl   FILE pointer of GRIB stream file which was opened by
!             a C program, GOPEN
!
!  OUTPUT:
!
!    fvar     Floating array to be written into GRIB file
!    ivar     Integer  array to be written into GRIB file
!
!    wdlen    Length of machine word (i.e. size of integer) in bytes
!
!    ipds     Integer array to carry the parameters for generating
!             PDS (Section 1).
!    igds     Integer array to carry the parameters for generating
!             GDS (Section 3).
!
!    nbufsz   Size of buffer of a GRIB message array
!    bds      BDS (GRIB Section 4)
!    mgrib    Working buffer of GRIB message array
!
!    nbtyes   Position at which the next GRIB message starts
!
!  WORK ARRAY:
!
!    tem1     Temporary work array.
!    item1    Integer temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in x,y,z dir.

  INTEGER :: wstag             ! Staggering flag for w-grid

  REAL :: z(nz)                ! Vertical coordinates

  INTEGER :: dflag             ! Dimension flag, 0: 3-D, 1: 2-D

  INTEGER :: nchanl            ! FILE pointer indicates the GRIB file

  INTEGER :: wdlen             ! Length of machine word

  INTEGER :: ipds(25)          ! ipds
  INTEGER :: igds(25)          ! igds
  INTEGER :: ibdshd(4)         ! BDS header

  INTEGER :: nbufsz            ! Size of GRIB message array
  CHARACTER (LEN=1) :: mgrib(nbufsz) ! GRIB message

  REAL :: fvar(nx,ny,nz)    ! Floating array
  INTEGER :: ivar(nx,ny,nz)    ! Integer  array

  CHARACTER (LEN=1) :: bds(nbufsz)   ! BDS
  INTEGER :: ibds(nbufsz/wdlen) ! Identical to BDS

  REAL :: tem1 (nx,ny)         ! Temporary work array
  INTEGER :: item1(nx,ny)      ! Integer temporary work array
!
!-----------------------------------------------------------------------
!
!  GRIB parameters
!
!-----------------------------------------------------------------------
!
  INTEGER :: msglen    ! Length of GRIB message
  INTEGER :: npts
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k, wstag1
  REAL :: tema
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  npts = nx*ny

  IF ( dflag == 1 ) THEN
    READ (nchanl) (mgrib(i),i=1,8)

    msglen = ICHAR(mgrib(5))*65536                                      &
           + ICHAR(mgrib(6))*256                                        &
           + ICHAR(mgrib(7))

    BACKSPACE (nchanl)

    READ (nchanl) (mgrib(i),i=1,msglen)

    CALL gribdec(npts,nz,                                               &
                 wdlen,ipds,igds,ibdshd,nbufsz,mgrib,                   &
                 tem1,item1, bds,ibds)

    DO j = 1,ny
      DO i = 1,nx
        fvar(i,j,1) = tem1 (i,j)
        ivar(i,j,1) = item1(i,j)
      END DO
    END DO

    RETURN
  END IF

  IF ( wstag > 1 .OR. wstag < 0 ) THEN
    WRITE(6,'(a,i1)')                                                   &
        'Do not know the w-grid type, wstag = ',wstag,                  &
        'Reset it to scalar grid.'
    wstag1 = 0
  ELSE
    wstag1 = wstag
  END IF

  DO k = 1,nz-1

    IF ( wstag1 == 0 ) THEN
      tema = 0.5 * (z(k) + z(k+1))
    ELSE
      tema = z(k)
    END IF

    ipds(11) = nint(tema)

    READ (nchanl) (mgrib(i),i=1,8)

    msglen = ICHAR(mgrib(5))*65536                                      &
           + ICHAR(mgrib(6))*256                                        &
           + ICHAR(mgrib(7))

    BACKSPACE (nchanl)

    READ (nchanl) (mgrib(i),i=1,msglen)

    CALL gribdec(npts,nz,                                               &
                 wdlen,ipds,igds,ibdshd,nbufsz,mgrib,                   &
                 tem1,item1, bds,ibds)

    DO j = 1,ny
      DO i = 1,nx
        fvar(i,j,k) = tem1 (i,j)
        ivar(i,j,k) = item1(i,j)
      END DO
    END DO

  END DO

  ipds(11) = nint(
  READ (nchanl) (mgrib(i),i=1,8)

  msglen = ICHAR(mgrib(5))*65536                                        &
         + ICHAR(mgrib(6))*256                                          &
         + ICHAR(mgrib(7))

  BACKSPACE (nchanl)

  READ (nchanl) (mgrib(i),i=1,msglen)

  CALL gribdec(npts,nz,                                                 &
               wdlen,ipds,igds,ibdshd,nbufsz,mgrib,                     &
               tem1,item1, bds,ibds)

  DO j = 1,ny
    DO i = 1,nx
      fvar(i,j,nz) = tem1 (i,j)
      ivar(i,j,nz) = item1(i,j)
    END DO
  END DO

  RETURN
END SUBROUTINE rdgrib
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE MKIPDS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE mkipds(nx,ny,nz, ipds) 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Make integer product data array ipds
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/27/1995
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!  OUTPUT:
!
!    ipds     IPDS array
!
!  WORK ARRAY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  INTEGER :: ipds(25)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ipds(1) = 28          ! number of bytes in PDS section
  ipds(2) = 0           ! parameter table version number, unknown
  ipds(3) = 0           ! ID number of originating center. CAPS?
  ipds(4) = 0           ! ID number of model. ARPS/CAPS?
  ipds(5) = 255         ! grid ID, 255 for self-defined GDS
  ipds(6) = 1           ! GDS flag, 1 = GDS section included
  ipds(7) = 0           ! BMS flag, 0 = no BMS section included
  ipds(8) = 0           ! indicator of parameter and unit
                        ! (table 2), depends on variables
  ipds(9) = 103         ! indicator of type of level (table 3)
                        ! depends on which variable
                        ! 103 for z coordinates in meters
                        ! 111 for soil layers in centimeters
! both use two octets: ipds(10) & ipds(11)
  ipds(10) = 0          ! value 1 of level, N/A
  ipds(11) = 0          ! value 2 of level (value of level)

  IF ( year <= 2000 ) THEN
    ipds(12) = year-1900 ! year of century
    ipds(23) = 20        ! century (20, change to 21 on Jan 1, 2001)
  ELSE
    ipds(12) = year-2000 ! year of century
    ipds(23) = 21        ! century (20, change to 21 on Jan 1, 2001)
  END IF

  ipds(13) = month      ! month of year
  ipds(14) = day        ! day of month
  ipds(15) = hour       ! hour of day
  ipds(16) = minute     ! minute of hour
  ipds(17) = 0          ! forecast time unit, minutes
  ipds(18) = 0          ! forecast time P1, or current time, curtim
  ipds(19) = 0          ! forecast time P2, N/A
  ipds(20) = 10         ! time range indicator,
                        ! 10 = use two octets from ipds(18)
  ipds(21) = 0          ! number include in average, N/A
  ipds(22) = 0          ! number missing from average, N/A
  ipds(24) = 0          ! subcenter ID number
  ipds(25) = 0          ! scaling power of 10,
                        ! depends on which variable

  RETURN
END SUBROUTINE mkipds
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE MKIGDS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE mkigds(nx,ny,nz,stagger,igds) 6,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Make integer grid description array IGDS.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/01/1995
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    stagger  Flag for staggered grid
!           = 0, s-grid, for scalar grid
!           = 1, u-grid, for u-vector grid
!           = 0, v-grid, for v-vector grid
!
!  OUTPUT:
!
!    igds     IPDS array
!
!  WORK ARRAY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  INTEGER :: stagger

  INTEGER :: igds(25)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i

  REAL :: swlat,swlon, nelat,nelon
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( stagger == 0 ) THEN
    swlat = swlats
    swlon = swlons
    nelat = nelats
    nelon = nelons
  ELSE IF ( stagger == 1 ) THEN
    swlat = swlatu
    swlon = swlonu
    nelat = nelatu
    nelon = nelonu
  ELSE IF ( stagger == 2 ) THEN
    swlat = swlatv
    swlon = swlonv
    nelat = nelatv
    nelon = nelonv
  ELSE
    WRITE(6,'(a)')                                                      &
        'Do not know the grid type. Set the grid type to scalar grid'
  END IF

  IF ( mapproj == 1 ) THEN       ! Polar stereographic
    igds(1) = nz        ! number of vertical (z) coordinates
    igds(2) = 255       ! location (in octets) of vertical
                        ! coordinate list
    igds(3) = 5         ! data representation type (table 6)
    igds(4) = nx        ! number of x-coordinates
    igds(5) = ny        ! number of y-coordinates
    igds(6) = nint(swlat*1000)  ! latitude  at southwest corner
    igds(7) = nint(swlon*1000)  ! longitude at southwest corner
    igds(8) = 8         ! u & v relative to x & y direction
    igds(9) = nint(trulon*1000) ! true longitude of projection
    igds(10) = dx
    igds(11) = dy

    IF ( trulat1 >= 0. ) THEN
      igds(12) = 0      ! northern pole on plane
    ELSE
      igds(12) = 1      ! southern pole on plane
    END IF

    igds(13) = 64       ! scanning flag, 64 for +i & +j direction
    DO i = 14,25
      igds(i) = 0       ! unused
    END DO

  ELSE IF ( mapproj == 2 .OR. mapproj == 0 ) THEN ! Lambert conformal
                                                  ! and no projection too

    igds( 1) = nz       ! number of vertical (z) coordinates
    igds( 2) = 255      ! location (in octets) of vertical
                        ! coordinate list
    igds( 3) = 3        ! data representation type (table 6)
    igds( 4) = nx       ! number of x-coordinates
    igds( 5) = ny       ! number of y-coordinates
    igds( 6) = nint(swlat*1000)  ! latitude  at southwest corner
    igds( 7) = nint(swlon*1000)  ! longitude at southwest corner
    igds( 8) = 8        ! u & v relative to x & y direction
    igds( 9) = nint(trulon*1000) ! true longitude of projection
    igds(10) = dx
    igds(11) = dy

    IF ( trulat1 >= 0. ) THEN
      igds(12) = 0      ! northern pole on plane
    ELSE
      igds(12) = 1      ! southern pole on plane
    END IF

    igds(13) = 64       ! scanning flag, 64 for +i & +j direction
    igds(14) = 0        ! unused
    igds(15) = nint(trulat1*1000)  ! first true latitude (millidegree)
    igds(16) = nint(trulat2*1000)  ! second true latitude (millidegree)
    igds(17) = -90000   ! latitude  of south pole (millidegree)
    igds(18) = nint(trulon*1000)   ! longitude of south pole (millidegree)

    DO i = 19, 25
      igds(i) = 0       ! unused
    END DO

  ELSE IF ( mapproj == 3 ) THEN  ! Mercator

    igds( 1) = nz       ! number of vertical (z) coordinates
    igds( 2) = 255      ! location (in octets) of vertical
                        ! coordinate list
    igds( 3) = 1        ! data representation type (table 6)
    igds( 4) = nx       ! number of x-coordinates
    igds( 5) = ny       ! number of y-coordinates
    igds( 6) = nint(swlat*1000)    ! latitude  at southwest corner
    igds( 7) = nint(swlon*1000)    ! longitude at southwest corner
    igds( 8) = 8        ! u & v relative to x & y direction
    igds( 9) = nint(nelat*1000)    ! latitude  of last point
    igds(10) = nint(nelon*1000)    ! longitude of last point
    igds(11) = nint((nelat-swlat)*1000./ny)  ! latitude  increment
    igds(12) = nint((nelon-swlon)*1000./nx)  ! longitude increment
    igds(13) = nint(trulat1*1000)  ! true latitude
    igds(14) = 64       ! scanning flag, 64 for +i & +j direction

    DO i = 15,25
      igds(i) = 0       ! unused
    END DO

  ELSE
    WRITE (6,'(a,i2/a)')                                                &
        'ARPS does support the type of map projection: ',mapproj,       &
        'Model stopped in MKIGDS.'
    CALL arpsstop(' ',1)
  END IF

  RETURN
END SUBROUTINE mkigds
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE GRIBCNTL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE gribcntl(nx,ny,nz,x,y,z,zp,gdsid, temxy1,temxy2) 1,6
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  ARPS GRIB data description file for GrADS display
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  11/27/1995
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    x        The x-coord. of the physical and computational grid.
!             Defined at u-point.
!    y        The y-coord. of the physical and computational grid.
!             Defined at v-point.
!    z        The z-coord. of the computational grid.
!             Defined at w-point on the staggered grid.
!
!  WORK ARRAY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  INTEGER :: gdsid             ! ID number of GRIB grid.
                               ! (should always be 255)

  REAL :: temxy1(nx,ny)        ! Temporary array
  REAL :: temxy2(nx,ny)        ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=120) :: grbctlfl
  CHARACTER (LEN=120) :: grbdatfl
  CHARACTER (LEN=15) :: chrstr, chr1
  CHARACTER (LEN=50) :: vartit(50), varnam(50)

  INTEGER :: varlev(50), varparam(50,4)

  INTEGER :: varnum
  INTEGER :: nchout0

  CHARACTER (LEN=3) :: monnam(12)            ! Name of months
  DATA monnam/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',                 &
              'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/

  CHARACTER (LEN=2) :: dtunit
  INTEGER :: ntm,tinc

  REAL :: lonmin,latmin,lonmax,latmax
  REAL :: xbgn,ybgn,zbgn
  REAL :: xinc,yinc,z0
  REAL :: lat11,lat12,lon11,lon12,lat21,lat22,lon21,lon22
  REAL :: latinc,loninc

  INTEGER :: fnlen
  INTEGER :: i,k, is
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Open the GrADS data control file: runname.ctl
!
!-----------------------------------------------------------------------
!
  IF ( mapproj /= 2 ) RETURN

  fnlen = lfnkey + 9
  grbctlfl(1:fnlen) = runname(1:lfnkey)//'.gribcntl'

  CALL fnversn( grbctlfl, fnlen )
  CALL getunit (nchout0)
  WRITE (6,'(a)') 'The GrADS data control file is '                     &
                    //grbctlfl(1:fnlen)

  OPEN (nchout0, FILE = grbctlfl(1:fnlen), STATUS = 'unknown')

  IF ( ldirnam == 0 ) THEN
    ldirnam = ldirnam + 1
    dirname(1:ldirnam) = '.'
  END IF

  fnlen = ldirnam + lfnkey + 14
  grbdatfl(1:fnlen) = dirname(1:ldirnam)//'/'                           &
                    //runname(1:lfnkey)//'.grb000000'

  xbgn = 0.5 * (x(1) + x(2))
  ybgn = 0.5 * (y(1) + y(2))
  zbgn = 0.5 * (z(1) + z(2))

  xinc = (x(2) - x(1))
  yinc = (y(2) - y(1))
  z0 = (z(2)-z(1))/2.

  CALL xytoll(nx,ny,x,y,temxy1,temxy2)

  CALL xytoll(1,1,xbgn,ybgn,lat11,lon11)

  CALL a3dmax0(temxy1,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1,                    &
               latmax,latmin)
  CALL a3dmax0(temxy2,1,nx,1,nx,1,ny,1,ny-1,1,1,1,1,                    &
               lonmax,lonmin)

  latinc = (latmax-latmin)/(ny-1)
  loninc = (lonmax-lonmin)/(nx-1)

  WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)')                                 &
           'latmin:latmax:latinc = ',                                   &
            latmin,':',latmax,':',latinc
  WRITE (6,'(a,f10.4,a,f10.4,a,f10.4)')                                 &
           'lonmin:lonmax:loninc = ',                                   &
           lonmin,':',lonmax,':',loninc

  WRITE (chrstr,'(i2.2,a,i2.2,a,i2.2,a3,i4.4)')                         &
      hour,':',minute,'Z',day,monnam(month),year

  IF ( nint(thisdmp) <= 0 ) THEN
    ntm = 1
  ELSE
    ntm = nint(tstop/thisdmp) + 1
  END IF

  IF ( nint(thisdmp) < 60 ) THEN
    WRITE (6, '(/a/a)')                                                 &
        'GrADS reqiures the smallest uint minute for time interval.',   &
        'Return to the caller.'
    tinc = nint(thisdmp)
    RETURN
  ELSE IF (thisdmp < 3600.) THEN
    tinc = nint(thisdmp/60.)
    dtunit = 'MN'
  ELSE IF (thisdmp < 86400.) THEN
    tinc = nint(thisdmp/3600.)
    dtunit = 'HR'
  ELSE
    tinc = nint(thisdmp/86400.)
    dtunit = 'DY'
  END IF

  varnum = 0

  IF ( varout == 1 ) THEN

    varnum = varnum + 1
    varnam(varnum) = 'u'
    vartit(varnum) = 'X-velocity total wind u (m/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 33
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'v'
    vartit(varnum) = 'Y-velocity total wind v (m/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 34
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'w'
    vartit(varnum) = 'Z-velocity total wind w (m/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 40
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'pt'
    vartit(varnum) = 'Potential Temperature (K)'
    varlev(varnum) = nz
    varparam(varnum,1) = 13
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'p'
    vartit(varnum) = 'Pressure (mb)'
    varlev(varnum) = nz
    varparam(varnum,1) = 1
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

  END IF

  IF ( mstout == 1 ) THEN

    varnum = varnum + 1
    varnam(varnum) = 'qv'
    vartit(varnum) = 'Water Vapor Mixing Ratio (g/kg)'
    varlev(varnum) = nz
    varparam(varnum,1) = 51
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'qc'
    vartit(varnum) = 'Cloud Water Mixing Ratio (g/kg)'
    varlev(varnum) = nz
    varparam(varnum,1) = 153
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'qr'
    vartit(varnum) = 'Rain Water Mixing Ratio (g/kg)'
    varlev(varnum) = nz
    varparam(varnum,1) = 170
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    IF ( iceout == 1 ) THEN

      varnum = varnum + 1
      varnam(varnum) = 'qi'
      vartit(varnum) = 'Cloud Ice Mixing Ratio (g/kg)'
      varlev(varnum) = nz
      varparam(varnum,1) = 178
      varparam(varnum,2) = 103
      varparam(varnum,3) = z0

      varnum = varnum + 1
      varnam(varnum) = 'qs'
      vartit(varnum) = 'Snow Mixing Ratio (g/kg)'
      varlev(varnum) = nz
      varparam(varnum,1) = 171
      varparam(varnum,2) = 103
      varparam(varnum,3) = z0

      varnum = varnum + 1
      varnam(varnum) = 'qh'
      vartit(varnum) = 'Hail Mixing Ratio (g/kg)'
      varlev(varnum) = nz
      varparam(varnum,1) = 179
      varparam(varnum,2) = 103
      varparam(varnum,3) = z0

    END IF

    IF ( rainout == 1 ) THEN

      varnum = varnum + 1
      varnam(varnum) = 'raing'
      vartit(varnum) = 'Grid Supersaturation Rain '
      varlev(varnum) = 0
      varparam(varnum,1) = 62
      varparam(varnum,2) = 1
      varparam(varnum,3) = 0

      varnum = varnum + 1
      varnam(varnum) = 'rainc'
      vartit(varnum) = 'Cumulus Convection Rain '
      varlev(varnum) = 0
      varparam(varnum,1) = 63
      varparam(varnum,2) = 1
      varparam(varnum,3) = 0

    END IF

    IF ( prcout == 1 ) THEN

      varnum = varnum + 1
      varnam(varnum) = 'prcrt1'
      vartit(varnum) = 'Total precipitation rate '
      varlev(varnum) = 1
      varparam(varnum,1) = 59
      varparam(varnum,2) = 111
      varparam(varnum,3) = 1

      varnum = varnum + 1
      varnam(varnum) = 'prcrt2'
      vartit(varnum) = 'Grid scale precipitation rate '
      varlev(varnum) = 2
      varparam(varnum,1) = 59
      varparam(varnum,2) = 111
      varparam(varnum,3) = 2

      varnum = varnum + 1
      varnam(varnum) = 'prcrt3'
      vartit(varnum) = 'cumulus precipitation rate '
      varlev(varnum) = 3
      varparam(varnum,1) = 59
      varparam(varnum,2) = 111
      varparam(varnum,3) = 3

      varnum = varnum + 1
      varnam(varnum) = 'prcrt4'
      vartit(varnum) = 'Microphysics precipitation rate '
      varlev(varnum) = 4
      varparam(varnum,1) = 59
      varparam(varnum,2) = 111
      varparam(varnum,3) = 4

    END IF
  END IF

  IF ( tkeout == 1 ) THEN

    varnum = varnum + 1
    varnam(varnum) = 'tke'
    vartit(varnum) = 'Turbulent kinetic energy (m**2/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 158
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

  END IF

  IF ( trbout == 1 ) THEN

    varnum = varnum + 1
    varnam(varnum) = 'kmh'
    vartit(varnum) = 'Horiz. Turb. Mixing Coeff. for '                  &
                     //'Momentum (m**2/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 185
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'kmv'
    vartit(varnum) = 'Vertical Turb. Mixing Coeff. for '                &
                     //'Momentum (m**2/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 186
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

  END IF

  IF ( sfcout == 1 ) THEN

    IF ( nstyp <= 1 ) THEN
      varnum = varnum + 1
      varnam(varnum) = 'ts'
      vartit(varnum) = 'Ground Surface Temperature (K)'
      varlev(varnum) = 0
      varparam(varnum,1) = 85
      varparam(varnum,2) = 111
      varparam(varnum,3) = 1

      varnum = varnum + 1
      varnam(varnum) = 't2'
      vartit(varnum) = 'Deep Soil Temperature (K)'
      varlev(varnum) = 0
      varparam(varnum,1) = 85
      varparam(varnum,2) = 111
      varparam(varnum,3) = 100

      varnum = varnum + 1
      varnam(varnum) = 'wg'
      vartit(varnum) = 'Surface Soil moisture'
      varlev(varnum) = 0
      varparam(varnum,1) = 144
      varparam(varnum,2) = 111
      varparam(varnum,3) = 1

      varnum = varnum + 1
      varnam(varnum) = 'w2'
      vartit(varnum) = 'Deep Soil moisture'
      varlev(varnum) = 0
      varparam(varnum,1) = 144
      varparam(varnum,2) = 111
      varparam(varnum,3) = 100

      varnum = varnum + 1
      varnam(varnum) = 'wr'
      vartit(varnum) = 'Canopy Soil moisture'
      varlev(varnum) = 0
      varparam(varnum,1) = 223
      varparam(varnum,2) = 111
      varparam(varnum,3) = 0

    ELSE

      DO is=0,nstyp
        WRITE (chr1,'(i1)') is

        varnum = varnum + 1
        varnam(varnum) = 'ts'//chr1
        vartit(varnum) = 'Ground Surface Temperature (K)'
        varlev(varnum) = 0
        varparam(varnum,1) = 85
        varparam(varnum,2) = 111
        varparam(varnum,3) = 1+is

        varnum = varnum + 1
        varnam(varnum) = 't2'//chr1
        vartit(varnum) = 'Deep Soil Temperature (K)'
        varlev(varnum) = 0
        varparam(varnum,1) = 85
        varparam(varnum,2) = 111
        varparam(varnum,3) = 100+is

        varnum = varnum + 1
        varnam(varnum) = 'wg'//chr1
        vartit(varnum) = 'Surface Soil moisture'
        varlev(varnum) = 0
        varparam(varnum,1) = 144
        varparam(varnum,2) = 111
        varparam(varnum,3) = 1+is

        varnum = varnum + 1
        varnam(varnum) = 'w2'//chr1
        vartit(varnum) = 'Deep Soil moisture'
        varlev(varnum) = 0
        varparam(varnum,1) = 144
        varparam(varnum,2) = 111
        varparam(varnum,3) = 100+is

        varnum = varnum + 1
        varnam(varnum) = 'wr'//chr1
        vartit(varnum) = 'Canopy Soil moisture'
        varlev(varnum) = 0
        varparam(varnum,1) = 223
        varparam(varnum,2) = 111
        varparam(varnum,3) = is
      END DO

    END IF

    IF ( snowout == 1 ) THEN
      varnum = varnum + 1
      varnam(varnum) = 'snowd'
      vartit(varnum) = 'Snow depth (m)'
      varlev(varnum) = 0
      varparam(varnum,1) = 66
      varparam(varnum,2) = 1
      varparam(varnum,3) = 0
    END IF

  END IF

  IF( radout == 1 ) THEN
    varnum = varnum + 1
    varnam(varnum) = 'radfrc'
    vartit(varnum) = 'Radiation forcing (K/s)'
    varlev(varnum) = nz
    varparam(varnum,1) = 216
    varparam(varnum,2) = 103
    varparam(varnum,3) = z0

    varnum = varnum + 1
    varnam(varnum) = 'radsw'
    vartit(varnum) = 'Net short wave radiation flux (W/m**2)'
    varlev(varnum) = 0
    varparam(varnum,1) = 204
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0

    varnum = varnum + 1
    varnam(varnum) = 'rnflx'
    vartit(varnum) = 'Total radiation flux (W/m**2)'
    varlev(varnum) = 0
    varparam(varnum,1) = 232
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0
  END IF

  IF( flxout == 1 ) THEN
    varnum = varnum + 1
    varnam(varnum) = 'uflx'
    vartit(varnum) = 'u flux'
    varlev(varnum) = 0
    varparam(varnum,1) = 124
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0

    varnum = varnum + 1
    varnam(varnum) = 'vflx'
    vartit(varnum) = 'v flux'
    varlev(varnum) = 0
    varparam(varnum,1) = 125
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0

    varnum = varnum + 1
    varnam(varnum) = 'ptflx'
    vartit(varnum) = 'pt flux'
    varlev(varnum) = 0
    varparam(varnum,1) = 122
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0

    varnum = varnum + 1
    varnam(varnum) = 'qvflx'
    vartit(varnum) = 'qv flux'
    varlev(varnum) = 0
    varparam(varnum,1) = 121
    varparam(varnum,2) = 1
    varparam(varnum,3) = 0
  END IF

  DO i=1,varnum
    varparam(i,4) = 0
  END DO

  WRITE (nchout0,'(a/a)')                                               &
      'TITLE   ARPS 4.0 GRIB Control for GrADS Display of '             &
      //runname(1:lfnkey)//' (a sample)','*'

  WRITE (nchout0,'(a,a)')                                               &
      'DSET    ',grbdatfl(1:fnlen)

  WRITE (nchout0,'(a,i4)')                                              &
      'DTYPE   GRIB ',gdsid

  WRITE (nchout0,'(a)')                                                 &
      'OPTIONS template'

  WRITE (nchout0,'(a)')                                                 &
      'INDEX   '//runname(1:lfnkey)//'.gribmap'

  WRITE (nchout0,'(a/a)')                                               &
      'UNDEF   -9.e+33','*'

  IF ( mapproj == 2 ) THEN

    WRITE (nchout0,'(a)')                                               &
        '* For lat-lon-lev display, umcomment the following 4 lines.'

    WRITE (nchout0,'(a,1x,i8,1x,i3,a,2f12.6,2i3,3f12.6,2f12.2)')        &
        'PDEF',nx,ny,' LCC',lat11,lon11,1,1,                            &
              trulat1,trulat2,trulon,xinc,yinc

  END IF

  WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)')                          &
      'XDEF',nx, '  LINEAR  ',lonmin,loninc

  WRITE (nchout0,'(a,1x,i8,a,f10.4,1x,f10.4)')                          &
      'YDEF',ny, '  LINEAR  ',latmin,latinc

  WRITE (nchout0,'(a,1x,i8,a)')                                         &
      'ZDEF',nz,'  LEVELS  '
  WRITE (nchout0,'(8f10.2)')                                            &
      ((z(k)+z(k+1))/2.,k=1,nz-1),z(nz)
  WRITE (nchout0,'(a/a/a)')                                             &
  '* WARNING& ! The vertical levels were set to computational ',        &
  '*          coordinates, z(k), because zp is not uniform in ',        &
  '*          horizontal when terrain option was turned on'

  WRITE (nchout0,'(a,1x,i8,a,a,3x,i2.2,a/a)')                           &
      'TDEF',ntm,'  LINEAR  ',chrstr,tinc,dtunit,'*'

  WRITE (nchout0,'(a,1x,i3)')                                           &
      'VARS',varnum

  DO i = 1, varnum
    WRITE (nchout0,'(a6,1x,i3,2x,i3.3,a,i4.4,a,i3.3,2x,a)')             &
        varnam(i),varlev(i), varparam(i,1),',',varparam(i,2),',',       &
        varparam(i,3),vartit(i)
  END DO

  WRITE (nchout0,'(a)')                                                 &
      'ENDVARS'

  CLOSE (nchout0)
  CALL retunit(nchout0)

  RETURN
END SUBROUTINE gribcntl