PROGRAM arps2eta212,139
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM ARSP2ETA212                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Converts ARPS history files to a ETA GRIB #212 format file.
!
!  Reads in a history file produced by ARPS in any ARPS format.
!  Converts variables and interpolates to a fixed set of pressure
!  levels.
!
!  NOTE:
!    Please make sure the ARPS domain covers ETA Grid #212 domain
!    before running this program.
!
!  INPUT Namelist
!    arps2eta212.input
!
!  Libraries:
!    libarps.a
!    libadas.a
!    HDF 4 library (When ARPS is in HDF 4 format)
!
!  Subroutine CALLs:
!    dtaread
!    mkarps2d
!    v2dint
!
!  Subroutine defined in this file
!    extrph
!    extrpt
!    extrpq
!    extrpuv
!    cal_avor
!    gtsinlat_ext
!    mkipds_212
!    mkigds_212
!    wrtgb        --  pack and write GRIB message for Grid #212
!    
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Yunheng Wang
!  04/19/2003
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  DATA ARRAYS READ IN:
!
!    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)
!    zpsoil   z coordinate of soil model (m) 
!
!    w        vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!    uprt     perturbation x velocity component (m/s)
!    vprt     perturbation y velocity component (m/s)
!    wprt     perturbation z velocity component (m/s)
!
!    qv       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)
!
!    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 density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsoil    Soil Temperature (K)
!    qsoil    Soil moisture 
!    wetcanp  Canopy water amount
!
!    rain     Total rain (raing + rainc)
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz          ! Grid dimensions.
  INTEGER :: nzsoil            ! Soil levels 
  INTEGER :: nstyps            ! Maximum number of soil types.

  INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,lenfil
  PARAMETER (nhisfile_max=200)
  CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Arrays to be read in:
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: x(:)       ! The x-coord. of the physical and
                                  ! computational grid.
                                  ! Defined at u-point.
  REAL, ALLOCATABLE :: y(:)       ! The y-coord. of the physical and
                                  ! computational grid.
                                  ! Defined at v-point.
  REAL, ALLOCATABLE :: z(:)       ! The z-coord. of the computational
                                  ! grid. Defined at w-point.

  REAL, ALLOCATABLE :: zp(:,:,:)  ! The height of the terrain.
  REAL, ALLOCATABLE :: zpsoil(:,:,:)  ! The height of the soil model.

  REAL, ALLOCATABLE :: uprt   (:,:,:) ! Perturbation u-velocity (m/s)
  REAL, ALLOCATABLE :: vprt   (:,:,:) ! Perturbation v-velocity (m/s)
  REAL, ALLOCATABLE :: wprt   (:,:,:) ! Perturbation w-velocity (m/s)
  REAL, ALLOCATABLE :: ptprt  (:,:,:) ! Perturbation potential temperature (K)
  REAL, ALLOCATABLE :: pprt   (:,:,:) ! Perturbation pressure (Pascal)
  REAL, ALLOCATABLE :: qvprt  (:,:,:) ! Perturbation water vapor specific
                                      ! humidity (kg/kg)
  REAL, ALLOCATABLE :: qc     (:,:,:) ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qr     (:,:,:) ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qi     (:,:,:) ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qs     (:,:,:) ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qh     (:,:,:) ! Hail mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: qv     (:,:,:) ! Water vapor specific humidity (kg/kg)

  REAL, ALLOCATABLE :: tke    (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: kmh    (:,:,:) ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: kmv    (:,:,:) ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: ubar   (:,:,:) ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vbar   (:,:,:) ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: wbar   (:,:,:) ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: ptbar  (:,:,:) ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: pbar   (:,:,:) ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3)
  REAL, ALLOCATABLE :: qvbar  (:,:,:) ! Base state water vapor specific

  INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
  REAL,    ALLOCATABLE :: stypfrct(:,:,:) ! Soil type
  INTEGER, ALLOCATABLE :: vegtyp(:,:)     ! Vegetation type
  REAL, ALLOCATABLE :: lai    (:,:)   ! Leaf Area Index
  REAL, ALLOCATABLE :: roufns (:,:)   ! Surface roughness
  REAL, ALLOCATABLE :: veg    (:,:)   ! Vegetation fraction

  REAL, ALLOCATABLE :: tsoil  (:,:,:,:) ! Soil Temperature (K)
  REAL, ALLOCATABLE :: qsoil  (:,:,:,:) ! Soil Moisture 
  REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
  REAL, ALLOCATABLE :: snowdpth(:,:)  ! Snow depth (m)

  REAL, ALLOCATABLE :: rain (:,:)     ! Total rainfall
  REAL, ALLOCATABLE :: raing(:,:)     ! Grid supersaturation rain
  REAL, ALLOCATABLE :: rainc(:,:)     ! Cumulus convective rain
  REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (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, ALLOCATABLE :: radfrc(:,:,:)  ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: radsw (:,:)    ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: rnflx (:,:)    ! Net radiation flux absorbed by surface
  REAL, ALLOCATABLE :: radswnet(:,:)  ! Net shortwave radiation
  REAL, ALLOCATABLE :: radlwin(:,:)   ! Incoming longwave radiation

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

  REAL, ALLOCATABLE :: e_mb   (:,:,:)    ! vapor pressure in mb
  REAL, ALLOCATABLE :: mix    (:,:,:)    ! total mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: esat_mb(:,:,:)    ! saturation vapor pressure in mb
  REAL, ALLOCATABLE :: rh     (:,:,:)    ! Relative humidity in %
  REAL, ALLOCATABLE :: t_dew  (:,:,:)    ! dewpoint temp. in degrees K

  REAL, ALLOCATABLE :: raing3hr(:,:)     ! 3 hour precipitation accumulation
  REAL, ALLOCATABLE :: rainc3hr(:,:)     ! 3 hour precipitation accumulation
  REAL, ALLOCATABLE :: rain3hr(:,:)      ! 3 hour precipitation accumulation
  REAL, ALLOCATABLE :: raingaccum(:,:)   ! precipitation accumulation 3 hr before
  REAL, ALLOCATABLE :: raincaccum(:,:)   ! precipitation accumulation 3 hr before
!
!-----------------------------------------------------------------------
!
!  Misc ARPS variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nchr, nchw, nchoutr
  REAL :: latnot(2), swx,swy,ctrx,ctry
  LOGICAL :: init
  REAL :: rlnp, soildp

  REAL, ALLOCATABLE :: xs(:)
  REAL, ALLOCATABLE :: ys(:)
  REAL, ALLOCATABLE :: zps(:,:,:)
  REAL, ALLOCATABLE :: zps_km(:,:,:)

  REAL :: gamma,ex1,ex2,rln700,p00,t00
  REAL, ALLOCATABLE :: p_mb(:,:,:)
  
  REAL, ALLOCATABLE :: mf2d(:,:)

  REAL, ALLOCATABLE :: avor(:,:,:)

  REAL, ALLOCATABLE :: uear(:,:,:)
  REAL, ALLOCATABLE :: vear(:,:,:)
!
!-----------------------------------------------------------------------
!
!  2-D stability index arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: lcl(:,:)   ! lifting condensation level
  REAL, ALLOCATABLE :: lfc(:,:)   ! level of free convection
  REAL, ALLOCATABLE :: el(:,:)    ! equilibrium level
  REAL, ALLOCATABLE :: twdf(:,:)  ! max. wet bulb pot. temp. difference
  REAL, ALLOCATABLE :: li(:,:)    ! lifted index
  REAL, ALLOCATABLE :: pbe(:,:)   ! CAPE
  REAL, ALLOCATABLE :: mbe(:,:)   ! Moist CAPE
  REAL, ALLOCATABLE :: nbe(:,:)   ! CIN
  REAL, ALLOCATABLE :: tcap(:,:)  ! Cap Strength

  REAL, ALLOCATABLE :: heli(:,:)  ! helicity
  REAL, ALLOCATABLE :: brn(:,:)   ! Bulk Richardson Number (Weisman and Klemp)
  REAL, ALLOCATABLE :: brnu(:,:)  ! Shear parameter of BRN, "U"
  REAL, ALLOCATABLE :: srlfl(:,:) ! storm-relative low-level flow (0-2km AGL)
  REAL, ALLOCATABLE :: srmfl(:,:) ! storm-relative mid-level flow (2-9km AGL)
  REAL, ALLOCATABLE :: shr37(:,:) ! 7km - 3km wind shear
  REAL, ALLOCATABLE :: ustrm(:,:) ! Estimated storm motion (modified Bob Johns)
  REAL, ALLOCATABLE :: vstrm(:,:) ! Estimated storm motion (modified Bob Johns)
  REAL, ALLOCATABLE :: blcon(:,:) ! Boundary layer convergence
  
  REAL, ALLOCATABLE :: wrk1(:),wrk2(:),wrk3(:),wrk4(:),wrk5(:),wrk6(:)
  REAL, ALLOCATABLE :: wrk7(:),wrk8(:),wrk9(:),wrk10(:),wrk11(:),wrk12(:)
  REAL, ALLOCATABLE :: tem2d1(:,:)
  REAL, ALLOCATABLE :: tem2d2(:,:)
  REAL, ALLOCATABLE :: tem2d3(:,:)
  REAL, ALLOCATABLE :: tem2d4(:,:)

!
!-----------------------------------------------------------------------
!
!  ETA output variables
!
!-----------------------------------------------------------------------
!
!  Constants

  INTEGER, PARAMETER :: nx_eta = 185, ny_eta = 129, nz_eta = 39
  REAL,    PARAMETER :: dx_eta = 40635.25, dy_eta = 40635.25  ! Meters

  INTEGER :: nprlvl
  INTEGER :: iprlvl(nz_eta)
  INTEGER, PARAMETER :: iprbgn = 1000, iprend = 50, iprinc = -25  ! mb

  INTEGER, PARAMETER :: iproj_eta = 3
  INTEGER, PARAMETER :: iproj_eia = 2
  REAL, PARAMETER    :: scale_eta = 1.0
  REAL, PARAMETER    :: latsw   =   12.190, lonsw   = -133.459,         &
                        latne   =   57.290, lonne   =  -49.385,         &
                        lattru_eta(2) = (/25.00, 25.00/),               &
                        lontru_eta =   -95.00  
  REAL :: x0_eta, y0_eta

  REAL, ALLOCATABLE :: x_eta(:), y_eta(:)
  REAL, ALLOCATABLE :: lat_eta(:,:), lon_eta(:,:)
  REAL, ALLOCATABLE :: umap(:,:)
  REAL, ALLOCATABLE :: vmap(:,:)

!-----------------------------------------------------------------------
!
!  Working arrays
!
!-----------------------------------------------------------------------

  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)
  REAL, ALLOCATABLE :: tem4(:,:,:)
  REAL, ALLOCATABLE :: sinlat(:,:)
  
  REAL, ALLOCATABLE :: tem2d_eta(:,:)    ! store data to be packed and written
  REAL, ALLOCATABLE :: tem2d2_eta(:,:)   ! store data to be packed and written

  REAL, ALLOCATABLE    :: wrkhori(:,:,:) ! store horizontal interpolations
  REAL, ALLOCATABLE    :: wrkhori2(:,:,:)! store horizontal interpolations
  REAL, ALLOCATABLE    :: p_eta(:,:,:)   ! store horizontally interpolated pressure
  REAL, ALLOCATABLE    :: lnp_eta(:,:,:) ! store -ln(p) at ETA grid

  INTEGER, ALLOCATABLE :: iloc(:,:)      ! x-index/y-index location of 
  INTEGER, ALLOCATABLE :: jloc(:,:)      ! ETA grid point in the ARPS array

  REAL, ALLOCATABLE    :: x2d(:,:)       ! x/y coordinate of ETA grid point
  REAL, ALLOCATABLE    :: y2d(:,:)       ! in ARPS coordinate 

  REAL, ALLOCATABLE :: dxfld(:), dyfld(:)
  REAL, ALLOCATABLE :: rdxfld(:), rdyfld(:)
  REAL, ALLOCATABLE :: slopey(:,:),alphay(:,:),betay(:,:)

!-----------------------------------------------------------------------
!
! Grib message variables
!
!-----------------------------------------------------------------------
  INTEGER :: wdlen             ! Length of machine word
  INTEGER :: ipds(25)          ! PDS integer array
  INTEGER :: igds(25)          ! GDS integer array

!-----------------------------------------------------------------------
!
!  Variables to determine the length of machine word, wdlen
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=8) :: ctema,ctemb
  INTEGER :: itema,itemb
  EQUIVALENCE ( itema,ctema )
  EQUIVALENCE ( itemb,ctemb )
  DATA ctema / '12345678' /

  REAL,    ALLOCATABLE :: mtem(:,:)
  INTEGER, ALLOCATABLE :: mitem1(:,:), mitem2(:,:)
  

!-----------------------------------------------------------------------
!
!  Namelists
!
!-----------------------------------------------------------------------
!
  INTEGER :: iorder  ! order of polynomial for interpolation (1, 2 or 3)
  NAMELIST /interpolation_options/ iorder

  INTEGER :: vvelout
  NAMELIST /output/ dirname,grbpkbit,readyfl, vvelout
                     ! these variables are declared in globcst.inc
!
!-----------------------------------------------------------------------
!
!  External functions
!
!-----------------------------------------------------------------------
!
  REAL :: wmr2td,oe,dpt, compute_density
  EXTERNAL wmr2td,oe,dpt, compute_density
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,klev,ifile,ireturn
  INTEGER :: istatus, fcttm
  REAL :: time
  INTEGER :: p1time,p2time
  CHARACTER(LEN=132) :: filname
  CHARACTER(LEN=132) :: filnamr
  INTEGER :: lenstr

  INTEGER :: varid
  REAL    :: rho

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Initializations
!
!-----------------------------------------------------------------------
!
!
  WRITE(6,'(9(/5x,a))')                                                 &
      '###############################################################',&
      '###############################################################',&
      '##                                                           ##',&
      '##                Welcome to ARPS2ETA212                     ##',&
      '## a program that reads in history files generated by ARPS   ##',&
      '## and produces a ETA GRIB #212 format file.                 ##',&
      '##                                                           ##',&
      '###############################################################',&
      '###############################################################'
!
!-----------------------------------------------------------------------
!
!  Get the names of the input data files.
!
!-----------------------------------------------------------------------
!
  CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)

  lengbf = len_trim(grdbasfn)

  CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf),                   &
       nx,ny,nz,nzsoil,nstyps, ireturn)

  IF( ireturn /= 0 ) THEN
    PRINT*,'Problem occured when trying to get dimensions from data.'
    PRINT*,'Program stopped.'
    STOP
  END IF

  WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz

  nstyp = nstyps
!
!-----------------------------------------------------------------------
!
!  Get output options and the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
   grbpkbit = 16
   dirname  = './'
   readyfl  = 0
   vvelout  = 2      ! dump omega by default

   READ(5,interpolation_options)
   READ(5,output)

!-----------------------------------------------------------------------
!
! Allocate ARPS arrays
!
!-----------------------------------------------------------------------

  ALLOCATE(x      (nx))
  ALLOCATE(y      (ny))
  ALLOCATE(z      (nz))
  ALLOCATE(zp     (nx,ny,nz))
  ALLOCATE(zpsoil (nx,ny,nzsoil))
  ALLOCATE(zps_km(nx,ny,nz))
  ALLOCATE(p_mb(nx,ny,nz))

  ALLOCATE(uprt   (nx,ny,nz))
  ALLOCATE(vprt   (nx,ny,nz))
  ALLOCATE(wprt   (nx,ny,nz))
  ALLOCATE(ptprt  (nx,ny,nz))
  ALLOCATE(pprt   (nx,ny,nz))
  ALLOCATE(qvprt  (nx,ny,nz))
  ALLOCATE(qv     (nx,ny,nz))
  ALLOCATE(qc     (nx,ny,nz))
  ALLOCATE(qr     (nx,ny,nz))
  ALLOCATE(qi     (nx,ny,nz))
  ALLOCATE(qs     (nx,ny,nz))
  ALLOCATE(qh     (nx,ny,nz))
  ALLOCATE(tke    (nx,ny,nz))
  ALLOCATE(kmh    (nx,ny,nz))
  ALLOCATE(kmv    (nx,ny,nz))
  ALLOCATE(ubar   (nx,ny,nz))
  ALLOCATE(vbar   (nx,ny,nz))
  ALLOCATE(wbar   (nx,ny,nz))
  ALLOCATE(ptbar  (nx,ny,nz))
  ALLOCATE(pbar   (nx,ny,nz))
  ALLOCATE(rhobar (nx,ny,nz))
  ALLOCATE(qvbar  (nx,ny,nz))

  ALLOCATE(soiltyp (nx,ny,nstyps))
  ALLOCATE(stypfrct(nx,ny,nstyps))
  ALLOCATE(vegtyp (nx,ny))
  ALLOCATE(lai    (nx,ny))
  ALLOCATE(roufns (nx,ny))
  ALLOCATE(veg    (nx,ny))

  ALLOCATE(tsoil  (nx,ny,nzsoil,0:nstyps))
  ALLOCATE(qsoil  (nx,ny,nzsoil,0:nstyps))
  ALLOCATE(wetcanp(nx,ny,0:nstyps))
  ALLOCATE(snowdpth(nx,ny))

  ALLOCATE(rain (nx,ny))
  ALLOCATE(raing(nx,ny))
  ALLOCATE(rainc(nx,ny))
  ALLOCATE(prcrate(nx,ny,4))

  ALLOCATE(raing3hr(nx,ny))
  ALLOCATE(rainc3hr(nx,ny))
  ALLOCATE(rain3hr(nx,ny))
  ALLOCATE(raingaccum(nx,ny))
  ALLOCATE(raincaccum(nx,ny))

  ALLOCATE(radfrc(nx,ny,nz))
  ALLOCATE(radsw (nx,ny))
  ALLOCATE(rnflx (nx,ny))
  ALLOCATE(radswnet(nx,ny))
  ALLOCATE(radlwin(nx,ny))

  ALLOCATE(usflx (nx,ny))
  ALLOCATE(vsflx (nx,ny))
  ALLOCATE(ptsflx(nx,ny))
  ALLOCATE(qvsflx(nx,ny))

  x      =0.0
  y      =0.0
  z      =0.0
  zp     =0.0

  uprt   =0.0
  vprt   =0.0
  wprt   =0.0
  ptprt  =0.0
  pprt   =0.0
  qvprt  =0.0
  qc     =0.0
  qr     =0.0
  qi     =0.0
  qs     =0.0
  qh     =0.0
  tke    =0.0
  kmh    =0.0
  kmv    =0.0
  ubar   =0.0
  vbar   =0.0
  wbar   =0.0
  ptbar  =0.0
  pbar   =0.0
  rhobar =0.0
  qvbar  =0.0
  qv     =0.0

  soiltyp =0.0
  stypfrct=0.0
  vegtyp =0.0
  lai    =0.0
  roufns =0.0
  veg    =0.0

  tsoil  =0.0
  qsoil  =0.0
  wetcanp=0.0
  snowdpth=0.0

  rain =0.0
  raing=0.0
  rainc=0.0
  prcrate=0.0

  raing3hr = 0.0
  rainc3hr = 0.0
  rain3hr  = 0.0
  raingaccum = 0.0
  raingaccum = 0.0

  radfrc=0.0
  radsw =0.0
  rnflx =0.0
  radswnet =0.0
  radlwin =0.0

  usflx =0.0
  vsflx =0.0
  ptsflx=0.0
  qvsflx=0.0

  ALLOCATE(xs(nx), STAT = istatus)
  ALLOCATE(ys(ny), STAT = istatus)
  ALLOCATE(zps(nx,ny,nz), STAT = istatus)
  zps    = 0.0
  zps_km = 0.0 
  p_mb   = 0.0

  ALLOCATE ( uear  (nx,ny,nz))     
  ALLOCATE ( vear  (nx,ny,nz))     
  uear = 0.0
  vear = 0.0

  ALLOCATE ( e_mb   (nx,ny,nz))    ! vapor pressure in mb

  ALLOCATE ( mix    (nx,ny,nz))    ! total mixing ratio (kg/kg)
  ALLOCATE ( esat_mb(nx,ny,nz))    ! saturation vapor pressure in mb
  ALLOCATE ( rh     (nx,ny,nz))    ! Relative humidity in %
  ALLOCATE ( t_dew  (nx,ny,nz))    ! dewpoint temp. in degrees K

  ALLOCATE ( avor  (nx,ny,nz))     ! absolute vorticity
  avor = 0.0
  
  ALLOCATE(lcl(nx,ny))
  ALLOCATE(lfc(nx,ny))
  ALLOCATE(el(nx,ny))
  ALLOCATE(twdf(nx,ny))
  ALLOCATE(li(nx,ny))
  ALLOCATE(pbe(nx,ny))
  ALLOCATE(mbe(nx,ny))
  ALLOCATE(nbe(nx,ny))
  ALLOCATE(tcap(nx,ny))
  lcl =0.0
  lfc =0.0
  el  =0.0
  twdf=0.0
  li  =0.0
  pbe =0.0
  mbe =0.0
  nbe =0.0
  tcap=0.0
  
  ALLOCATE(wrk1(nz),wrk2(nz),wrk3(nz),wrk4(nz),wrk5(nz),wrk6(nz))
  ALLOCATE(wrk7(nz),wrk8(nz),wrk9(nz),wrk10(nz),wrk11(nz),wrk12(nz))
  
  wrk1 = 0.0
  wrk2 = 0.0
  wrk3 = 0.0
  wrk4 = 0.0
  wrk5 = 0.0
  wrk6 = 0.0
  wrk7 = 0.0
  wrk8 = 0.0
  wrk9 = 0.0
  wrk10= 0.0
  wrk11= 0.0
  wrk12= 0.0
  
  ALLOCATE(tem2d1(nx,ny))
  ALLOCATE(tem2d2(nx,ny))
  ALLOCATE(tem2d3(nx,ny))
  ALLOCATE(tem2d4(nx,ny))  
  tem2d1=0.0
  tem2d2=0.0
  tem2d3=0.0
  tem2d4=0.0  
  
  ALLOCATE( mf2d(nx,ny))
  mf2d=0.0
  
  ALLOCATE(heli (nx,ny))
  ALLOCATE(brn  (nx,ny))
  ALLOCATE(brnu (nx,ny))
  ALLOCATE(srlfl(nx,ny))
  ALLOCATE(srmfl(nx,ny))
  ALLOCATE(shr37(nx,ny))
  ALLOCATE(ustrm(nx,ny))
  ALLOCATE(vstrm(nx,ny))
  ALLOCATE(blcon(nx,ny))
    
  heli =0.0
  brn  =0.0
  brnu =0.0
  srlfl=0.0
  srmfl=0.0
  shr37=0.0
  ustrm=0.0
  vstrm=0.0
  blcon=0.0
!---------------------------------------------------------------------------
!
! Allocate eta arrays
!
!---------------------------------------------------------------------------

  ALLOCATE(x_eta  (nx_eta), STAT = istatus)  
  ALLOCATE(y_eta  (ny_eta), STAT = istatus)  
  ALLOCATE(lat_eta(nx_eta,ny_eta), STAT = istatus)  
  ALLOCATE(lon_eta(nx_eta,ny_eta), STAT = istatus)  

  ALLOCATE ( umap  (nx_eta,ny_eta))     
  ALLOCATE ( vmap  (nx_eta,ny_eta))     
  umap = 0.0
  vmap = 0.0

!--------------------------------------------------------------------------
!
! Allocate temporary arrays
!
!--------------------------------------------------------------------------

  ALLOCATE(lnp_eta(nx_eta,ny_eta,nz), STAT = istatus)
  lnp_eta = 0.0

  ALLOCATE(p_eta(nx_eta,ny_eta,nz), STAT = istatus)
  p_eta = 0.0

  ALLOCATE(wrkhori(nx_eta,ny_eta,nz), STAT = istatus)  
  wrkhori = 0.0
  ALLOCATE(wrkhori2(nx_eta,ny_eta,nz), STAT = istatus)  
  wrkhori2 = 0.0

  ALLOCATE(iloc(nx_eta,ny_eta), STAT = istatus)  
  ALLOCATE(jloc(nx_eta,ny_eta), STAT = istatus)  
  iloc = 0
  jloc = 0
  ALLOCATE(x2d (nx_eta,ny_eta), STAT = istatus)  
  ALLOCATE(y2d (nx_eta,ny_eta), STAT = istatus)  
  x2d = 0.0
  y2d = 0.0

  ALLOCATE( dxfld(nx), STAT = istatus)
  ALLOCATE( dyfld(ny), STAT = istatus)
  ALLOCATE(rdxfld(nx), STAT = istatus)
  ALLOCATE(rdyfld(ny), STAT = istatus)
  dxfld = 0.0
  dyfld = 0.0
  rdxfld = 0.0
  rdyfld = 0.0

  ALLOCATE(slopey(nx,ny), STAT = istatus)
  ALLOCATE(alphay(nx,ny), STAT = istatus)
  ALLOCATE( betay(nx,ny), STAT = istatus)
  slopey = 0.0
  alphay = 0.0
  betay  = 0.0

  ALLOCATE(tem1(nx,ny,nz))
  ALLOCATE(tem2(nx,ny,nz))
  ALLOCATE(tem3(nx,ny,nz))
  ALLOCATE(tem4(nx,ny,nz))
  ALLOCATE(sinlat(nx,ny))
  sinlat = 0.0
  
  tem1 =0.0
  tem2 =0.0
  tem3 =0.0
  tem4 =0.0

  ALLOCATE(tem2d_eta(nx_eta,ny_eta), STAT = istatus)
  tem2d_eta = 0.0
  ALLOCATE(tem2d2_eta(nx_eta,ny_eta), STAT = istatus)
  tem2d2_eta = 0.0

!-----------------------------------------------------------------------
!
!  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
    WRITE (6,'(a,i2,a)')  &
        'The length of machine word is ',wdlen,' bytes'

!-----------------------------------------------------------------------
!
! Initilize GRIB message variables
!
!-----------------------------------------------------------------------

  CALL mkigds_212(nx_eta,ny_eta,nprlvl,igds)

!
  ALLOCATE(mtem(nx_eta,ny_eta),   STAT = istatus)
  ALLOCATE(mitem1(nx_eta,ny_eta), STAT = istatus)
  ALLOCATE(mitem2(nx_eta,ny_eta), STAT = istatus)
  mtem = 0.0
  mitem1 = 0
  mitem2 = 0

!-----------------------------------------------------------------------
!
!  Find x,y locations of ETA grid.
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj_eia,scale_eta,lattru_eta,lontru_eta)
  CALL setorig(2,latsw,lonsw)
  CALL lltoxy(1,1,latsw,lonsw,x0_eta,y0_eta) 
!  CALL setorig(1,x0_eta,y0_eta)

  DO i=1,nx_eta
    x_eta(i)=x0_eta+(i-1)*dx_eta
  END DO
  DO j=1,ny_eta
    y_eta(j)=y0_eta+(j-1)*dy_eta
  END DO
  CALL xytoll(nx_eta,ny_eta,x_eta,y_eta,lat_eta,lon_eta)

  DO klev = iprbgn,iprend,iprinc
    iprlvl((iprbgn-klev)/ABS(iprinc) + 1) = klev
  END DO
  nprlvl = SIZE(iprlvl)
  IF(nprlvl > nz_eta) THEN
    WRITE(6,*) 'ERROR: nz_eta is too small.'
    STOP
  END IF

!--------------------------------------------------------------------
!
! Begin loop over all of the history files
!
!--------------------------------------------------------------------

  lengbf=LEN_trim(grdbasfn)
  WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)

  init = .false.
  nchw = 10

  DO ifile = 1,nhisfile
    lenfil = LEN(hisfile(ifile))
    CALL strlnth( hisfile(ifile), lenfil)
    WRITE(6,'(/a,/1x,a)')' The data set name is ',                      &
                                           hisfile(ifile)(1:lenfil)

!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
    CALL dtaread(nx,ny,nz,nzsoil,nstyps,hinfmt,nchr,grdbasfn(1:lengbf), &
               lengbf,                                                  &
               hisfile(ifile)(1:lenfil),lenfil,time,                    &
               x,y,z,zp,zpsoil,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,                  &
               tsoil,qsoil,wetcanp,snowdpth,                            &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,radswnet,radlwin,                     &
               usflx,vsflx,ptsflx,qvsflx,                               &
               ireturn, tem1,tem2,tem3)

!    IF(nx < nx_eta .OR. ny < ny_eta) THEN
!      WRITE(6,'(a)') 'ARPS domain is smaller than ETA 212 domain.'
!      STOP
!    END IF
!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
    IF( ireturn == 0 ) THEN   ! successful read

    !WRITE(6,*) 'time =', time

    fcttm = INT(time/3600)
    CALL gtlfnkey(runname, lenstr)

    WRITE(filname,'(a,a1,I4.4,I2.2,I2.2,I2.2,a,I2.2,a4)')              &
           runname(1:lenstr),'.',year,month,day,hour,'f',fcttm,'.grb'

    p1time = INT((time-10800)/3600)      ! in minutes
    p2time = INT(time/3600)              ! in minutes

    lenstr = LEN_TRIM(dirname)
    IF(lenstr > 0) THEN
      IF(dirname(lenstr:lenstr) /= '/') THEN
        dirname(lenstr+1:lenstr+1) = '/'
        lenstr = lenstr + 1
      END IF
      filname = dirname(1:lenstr) // TRIM(filname)
      print *, 'Output file name is ',filname
    END IF

!    lenstr = 132
!    CALL strlnth( filname, lenstr )
    lenstr = LEN_TRIM(filname)

    CALL asnctl ('NEWLOCAL', 1, istatus)
    CALL asnfile(filname, '-F f77 -N ieee', istatus)

    CALL getunit( nchw )

    OPEN(UNIT=nchw,FILE=trim(filname),STATUS='unknown',             &
                   FORM='unformatted',IOSTAT= istatus )

    IF(istatus /= 0) THEN
      WRITE(6,*) 'Error while creating file.',TRIM(filname)
      STOP
    END IF

    IF( .NOT.init ) THEN

!-----------------------------------------------------------------------
!
!  Establish coordinate for scalar fields.
!
!-----------------------------------------------------------------------
!
        DO i=1,nx-1
          xs(i)=0.5*(x(i)+x(i+1))
        END DO
        xs(nx)=2.*xs(nx-1)-xs(nx-2)

        DO j=1,ny-1
          ys(j)=0.5*(y(j)+y(j+1))
        END DO
        ys(ny)=2.*ys(ny-1)-ys(ny-2)

        DO k=1,nz-1
          DO j=1,ny-1
            DO i=1,nx-1
              zps(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
            END DO
          END DO
        END DO
!
!-----------------------------------------------------------------------
!
!  Find x,y locations of ETA grid in terms of the ARPS grid.
!
!-----------------------------------------------------------------------
        latnot(1) = trulat1
        latnot(2) = trulat2
        CALL setmapr(mapproj,1.0/sclfct,latnot,trulon)
        CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
        swx = ctrx - (FLOAT((nx-3))/2.) * dx/sclfct
        swy = ctry - (FLOAT((ny-3))/2.) * dy/sclfct
        CALL setorig(1,swx,swy)
        
        CALL xytomf(nx,ny,xs,ys,mf2d)
        
        CALL lltoxy(nx_eta,ny_eta,lat_eta,lon_eta,x2d,y2d)

        CALL setijloc(nx_eta,ny_eta,nx,ny,x2d,y2d,                    &
                    xs,ys,iloc,jloc)


        CALL setdxdy(nx,ny,1,nx,1,ny, xs,ys,dxfld,dyfld,rdxfld,rdyfld)

        init=.true.
     END IF  ! NOT init

! Restore to use ETA map projection

     CALL setmapr(iproj_eia,scale_eta,lattru_eta,lontru_eta)
     CALL setorig(2,latsw,lonsw)
!
!-----------------------------------------------------------------------
!
!    Begin ARPS data conversions
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            pprt(i,j,k)=pprt(i,j,k)+pbar(i,j,k)
            ptprt(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
            qvprt(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
            tem1(i,j,k)=0.5*(uprt(i,j,k)+ubar(i,j,k)+                 &
                        uprt(i+1,j,k)+ubar(i+1,j,k))
            tem2(i,j,k)=0.5*(vprt(i,j,k)+vbar(i,j,k)+                 &
                        vprt(i,j+1,k)+vbar(i,j+1,k))
            tem3(i,j,k)=0.5*(wprt(i,j,k)+wbar(i,j,k)+                 &
                        wprt(i,j,k+1)+wbar(i,j,k+1))
            qi(i,j,k)=qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
          END DO
        END DO
      END DO

!
!-----------------------------------------------------------------------
!
!  Swap wind data back into wind arrays
!
!-----------------------------------------------------------------------
!

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            uprt(i,j,k)=tem1(i,j,k)
            vprt(i,j,k)=tem2(i,j,k)
            wprt(i,j,k)=tem3(i,j,k)
          END DO
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!  Put temperature into tem2
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem2(i,j,k)=ptprt(i,j,k)*(pprt(i,j,k)/p0)**rddcp
          END DO
        END DO
      END DO

!-----------------------------------------------------------------------
!
!  Processing Pressure data
!
!-----------------------------------------------------------------------

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,pprt(:,:,k),p_eta(:,:,k),              &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
        lnp_eta(:,:,k) = -ALOG(p_eta(:,:,k))
      END DO

!-----------------------------------------------------------------------
!
!  Processing Temperature data
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,tem2(:,:,k),wrkhori(:,:,k),            &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
        ! here use mkarps2d in the reverse direction
        ! interpolate from ARPS grid to ETA grid instead of the other
      END DO

      PRINT *, ' Storing Temperature GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels
        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpt(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message
        CALL mkipds_212(11,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO

!  Find temperature in 2m
      DO k = 1, nz
        DO j = 1,ny_eta
          DO i = 1, nx_eta
            wrkhori2(i,j,k) = z(k)
          END DO
        END DO 
      END DO
      
      PRINT *, ' Storing Temperature at 2m.'

      CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                   wrkhori,wrkhori2,2,tem2d_eta)
      CALL mkipds_212(11,105,2,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)

!-----------------------------------------------------------------------
!
!  Process Specific Humidity (51)
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,          &
                   xs,ys,x2d,y2d,qvprt(:,:,k),wrkhori(:,:,k),         &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,  &
                   tem4)
      END DO

      PRINT *, ' Storing Specific Humidity GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,      &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpq(nx_eta,ny_eta,nz,wrkhori,p_eta,                  &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message

        CALL mkipds_212(51,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,      &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO
 
!-----------------------------------------------------------------------
!
!  Process U wind (33) and V wind (34)
!
!  Wind needs to be rotated to ETA grid
!
!-----------------------------------------------------------------------

! First rotate the wind to true-earth-relative in ARPS map projection
      CALL setmapr(mapproj,1.0/sclfct,latnot,trulon)
      CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
      swx = ctrx - (FLOAT((nx-3))/2.) * dx/sclfct
      swy = ctry - (FLOAT((ny-3))/2.) * dy/sclfct
      CALL setorig(1,swx,swy)
      CALL xytoll(nx,ny,xs,ys,tem2d1,tem2d2)
      DO k = 1,nz
        CALL uvmptoe(nx,ny,uprt(:,:,k),vprt(:,:,k),                 &
                     tem2d2,uear(:,:,k),vear(:,:,k))
      END DO

! Second rotate the true earth wind to ETA map projection
 
      CALL setmapr(iproj_eia,scale_eta,lattru_eta,lontru_eta)
      CALL setorig(2,latsw,lonsw)

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1

        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,        &
                   xs,ys,x2d,y2d,uear(:,:,k),wrkhori(:,:,k),        &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,&
                   tem4)

        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,        &
                   xs,ys,x2d,y2d,vear(:,:,k),wrkhori2(:,:,k),       &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,&
                   tem4)
      END DO

!-----------------------------------------------------------------------
!
!  10m u and v
!  Dan recommended to use k = 2 is 10m
!
!-----------------------------------------------------------------------

      CALL uvetomp(nx_eta,ny_eta,wrkhori(:,:,2),wrkhori2(:,:,2),    &
                   lon_eta,umap,vmap)

      PRINT *, ' Storing horizontal wind at 10m.'

      CALL mkipds_212(33,105,10,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,umap,nchw,wdlen,grbpkbit,            &
                   ipds,igds, mtem,mitem1,mitem2)

      CALL mkipds_212(34,105,10,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,vmap,nchw,wdlen,grbpkbit,            &
                   ipds,igds, mtem,mitem1,mitem2)

!------------------------------------------------------------------------

      PRINT *, ' Storing Horizontal wind GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)
        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori2,lnp_eta,rlnp,tem2d2_eta)

        CALL extrpuv(nx_eta,ny_eta,nz,wrkhori,wrkhori2,p_eta,           &
                     iprlvl(klev),tem2d_eta,tem2d2_eta)

        ! 3. Pack and write GRIB message

        CALL uvetomp(nx_eta,ny_eta,tem2d_eta,tem2d2_eta,lon_eta,umap,vmap)

        CALL mkipds_212(33,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,umap,nchw,wdlen,grbpkbit,              &
                   ipds,igds, mtem,mitem1,mitem2)

        CALL mkipds_212(34,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,vmap,nchw,wdlen,grbpkbit,              &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO

!-----------------------------------------------------------------------
!
!  Process Geopotential Heigth (gpm) (7)
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,zps(:,:,k),wrkhori(:,:,k),             &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      PRINT *, ' Storing Geopotential Height GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrph(nx_eta,ny_eta,nz,wrkhori,wrkhori,p_eta,             &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message

        CALL mkipds_212(07,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,        &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO

!-----------------------------------------------------------------------
!
!  Process W wind (40)   (ETA uses Pa S**-1 -- 39)
!
!-----------------------------------------------------------------------

      IF(vvelout == 2) THEN   ! convert W to Omega
        DO k=1,nz
          DO j=1,ny
            DO i=1,nx
              rho = compute_density(tem2(i,j,k),pprt(i,j,k))
              wprt(i,j,k)= -rho*g*wprt(i,j,k) 
            END DO
          END DO
        END DO
        varid = 39
      ELSE IF(vvelout == 1) THEN
        varid = 40
      ELSE
        varid = 0
      END IF

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      IF(varid /= 0) THEN

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,wprt(:,:,k),wrkhori(:,:,k),            &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      PRINT *, ' Storing W-wind GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpq(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message

        CALL mkipds_212(varid,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO

      END IF   ! varid /= 0


!-----------------------------------------------------------------------
!
!  Process soil temperature (85)
!    level type = 111   Depth below land surface
!    ETA use 112 layer between two depths below land surface
!
!-----------------------------------------------------------------------

      PRINT *, ' Storing soil temperature at layer below land surface.'

      DO k = 1,nzsoil
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,tsoil(:,:,k,0),tem2d_eta,              &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
        soildp = 100*(zpsoil(2,2,1)-zpsoil(2,2,k))   

        CALL mkipds_212(85,111,soildp,0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

      END DO

!-----------------------------------------------------------------------
!
!  Process soil moisture content (144)
!    level type = 111   Depth below land surface
!    ETA use 112 layer between two depths below land surface
!
!-----------------------------------------------------------------------

      PRINT *, ' Storing soil moisture GRIB data.'

      DO k = 1,nzsoil
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,qsoil(:,:,k,0),tem2d_eta,              &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
        soildp = 100*(zpsoil(2,2,1)-zpsoil(2,2,k))   
                              !   layer below land surface, cm
        CALL mkipds_212(144,111,soildp,0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

      END DO

!-----------------------------------------------------------------------
!
!  Water equiv. of accum snow depth (kg/m**2) -- 65
!
!----------------------------------------------------------------------

      PRINT *, ' Storing accum. snow depth.'

      tem1(:,:,1) = 100* snowdpth(:,:)

      CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                  xs,ys,x2d,y2d,tem1(:,:,1),tem2d_eta,                &
                  dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                  tem4)

      CALL mkipds_212(65,01,00,0,p2time,0,10,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

!----------------------------------------------------------------------
!
!  Convective precip (63)
!
!  rainc unit is mm
!  63 require kg m**-2  = rhowater (1000 kg m**-3) * rainc (mm) * 1.0E-3
!
!----------------------------------------------------------------------
      IF (time >= 10800 .AND. MOD(INT(time),10800) == 0) THEN
                       ! precipitation dump every 3 hours  
        raing3hr = raing - raingaccum
        rainc3hr = rainc - raincaccum
        raingaccum = raing
        raingaccum = rainc
        rain3hr = raing3hr + rainc3hr

        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                    xs,ys,x2d,y2d,rain3hr,tem2d_eta,                    &
                    dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                    tem4)
        print *, ' Storing total precipitation rainfall.'
        CALL mkipds_212(61,01,00,4,p1time,p2time,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                    xs,ys,x2d,y2d,raing3hr,tem2d_eta,                   &
                    dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                    tem4)
        print *, ' Storing large scale precipitation.'
        CALL mkipds_212(62,01,00,4,p1time,p2time,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                    xs,ys,x2d,y2d,rainc3hr,tem2d_eta,                   &
                    dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                    tem4)
        print *, ' Storing convective rainfall.'
        CALL mkipds_212(63,01,00,4,p1time,p2time,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

      END IF
!-----------------------------------------------------------------------
!
!  Output near-sfc data.
!  Simularity theory or something could be applied here
!  to make these data be valid at sfc instrument height,
!  for now we output level 2.
!
!-----------------------------------------------------------------------

      CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                  xs,ys,x2d,y2d,pprt(:,:,2),tem2d_eta,                &
                  dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                  tem4)
      print *, ' Storing near-sfc presure.'

      CALL mkipds_212(01,01,0,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)

!-----------------------------------------------------------------------
!
!  Geopotential Height (gpm)
!
!-----------------------------------------------------------------------

      CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                  xs,ys,x2d,y2d,zps(:,:,2),tem2d_eta,                 &
                  dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                  tem4)

      print *, ' Storing near-sfc Geopotential Height'

      CALL mkipds_212(07,01,00,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,       &
                 ipds,igds, mtem,mitem1,mitem2)

!-----------------------------------------------------------------------
!
! Surface Temperature (11)
!
!-----------------------------------------------------------------------

      CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                  xs,ys,x2d,y2d,tem2(:,:,2),tem2d_eta,                &
                  dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,   &
                  tem4)

      print *, ' Storing near-sfc temperature'

      CALL mkipds_212(11,01,00,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)

!----------------------------------------------------------------------
!
!  Calculation t_dew, rh etc.
!
!-----------------------------------------------------------------------

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            zps_km(i,j,k)=zps(i,j,k)/1000.0
            p_mb(i,j,k)=0.01*pprt(i,j,k)
            mix(i,j,k)=1000.0*(qvprt(i,j,k)/(1.-qvprt(i,j,k)))
            e_mb(i,j,k)=qvprt(i,j,k)*p_mb(i,j,k)/(qvprt(i,j,k)+         &
                        287.0/461.0)
            esat_mb(i,j,k)=6.1078*EXP((2.500780E6/461.0)*((1.0/273.15)- &
                          (1.0/tem2(i,j,k))))
            t_dew(i,j,k)=wmr2td(p_mb(i,j,k),mix(i,j,k))                 &
                       + 273.15
            rh(i,j,k)=100.0*(e_mb(i,j,k)/esat_mb(i,j,k))
            rh(i,j,k)=MAX(0.,rh(i,j,k))
          END DO
        END DO
      END DO

!-----------------------------------------------------------------------
!
!  Process Relative Humidity (52)
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,rh(:,:,k),wrkhori(:,:,k),              &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      PRINT *, ' Storing Relative Humidity GRIB data.'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpq(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message

        DO j = 1, ny_eta
           DO i = 1, nx_eta
             tem2d_eta(i,j) = MAX( MIN(tem2d_eta(i,j), 100.0), 0.0)
           END DO
        END DO

        CALL mkipds_212(52,100,iprlvl(klev),0,p2time,0,0,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

      END DO
 
!-----------------------------------------------------------------------
!
!  Process Dew-point temperature (17) at 850mb only
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,t_dew(:,:,k),wrkhori(:,:,k),           &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      rlnp = -ALOG(100.*FLOAT(850))

      PRINT *, ' Storing Dew-point temp. at pressure = ',850, 'mb.'

      ! 2. Vertically interpolate to pressure levels

      CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

      CALL extrpt(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                  850,tem2d_eta)

      ! 3. Pack and write GRIB message

      CALL mkipds_212(17,100,850,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)
 
!----------------------------------------------------------------------
!
!  Find dew-point temperature at 2m
!
!---------------------------------------------------------------------

      DO k = 1, nz
        DO j = 1,ny_eta
          DO i = 1, nx_eta
            wrkhori2(i,j,k) = z(k)
          END DO
        END DO 
      END DO
      
      CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,           &
                   wrkhori,wrkhori2,2,tem2d_eta)

      PRINT *, ' Storing dew-point temperature at 2m'

      CALL mkipds_212(17,105,2,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)

!-----------------------------------------------------------------------
!
!  Process Clould ice, (58) kg m**-2
!
!-----------------------------------------------------------------------

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        DO j = 1, ny-1
          DO i = 1, nx-1
            tem1(i,j,k) = rhobar(i,j,k) * qi(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)

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,tem1(:,:,k),wrkhori(:,:,k),            &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      PRINT *, ' Storing Could ice GRIB data'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))

        ! 2. Vertically interpolate to pressure levels

        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpq(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message

        CALL mkipds_212(58,100,iprlvl(klev),0,p2time,0,10,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)
      END DO
 
!-----------------------------------------------------------------------
!
!  Calculate abosolute vorticity
!
!-----------------------------------------------------------------------

      CALL cal_avor(avor,uprt,vprt,x,y,nx,ny,nz,1,0,omega,              &
                     sinlat,xs,ys,tem2d1)

      avor = avor * 1.0E-5

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      DO k = 1,nz-1
        CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,avor(:,:,k),wrkhori(:,:,k),            &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)
      END DO

      PRINT *, ' Storing absolut vor. GRIB data'

      DO klev=1,nprlvl
        rlnp = -ALOG(100.*FLOAT(iprlvl(klev)))
        
        ! 2. Vertically interpolate to pressure levels
        CALL v2dinta(nx_eta,ny_eta,nz,1,nx_eta,1,ny_eta,1,nz-1,         &
                      wrkhori,lnp_eta,rlnp,tem2d_eta)

        CALL extrpq(nx_eta,ny_eta,nz,wrkhori,p_eta,                     &
                    iprlvl(klev),tem2d_eta)

        ! 3. Pack and write GRIB message


        CALL mkipds_212(41,100,iprlvl(klev),0,p2time,0,10,ipds)
        CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,         &
                   ipds,igds, mtem,mitem1,mitem2)

        IF(istatus /= 0) THEN
           WRITE(6,*) 'Error when presure =', iprlvl(klev)
           STOP
        END IF

      END DO

!
!-----------------------------------------------------------------------
!
!  Put temperature into tem2
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem2(i,j,k)=ptprt(i,j,k)*                                   &
                          (pprt(i,j,k)/p0)**rddcp
            tem3(i,j,k)=-ALOG(pprt(i,j,k))
          END DO
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!    Calculate temperature (K) at ARPS grid points
!   and 700mb pressure level
!
!-----------------------------------------------------------------------
!
      rln700=-ALOG(70000.0)
      CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                       &
                      tem2,tem3,rln700,tem1)

!-----------------------------------------------------------------------
!
!    Calculate sea level pressure (mb)
!    Reduction method: Benjamin and Miller: 1990, MWR, vol.118, No.10,
!                   Page: 2100-2101
!
!-----------------------------------------------------------------------
!
      gamma=.0065             ! std lapse rate per meter
      ex1=0.1903643           ! R*gamma/g
      ex2=5.2558774
      DO j=1,ny-1
        DO i=1,nx-1
          p00 = 0.01*(pprt(i,j,2))
          IF (p00 <= 700.0) THEN
            t00 = tem2(i,j,2)
          ELSE
            t00 = tem1(i,j,1)*(p00/700.0)**ex1
          END IF
          tem1(i,j,1)=p00*(1.0+gamma*zps(i,j,2)/t00)**ex2
        END DO
      END DO

      PRINT *, ' Storing MSL Pressure '

      ! 1. Horizontally interpolate from ARPS grid to ETA grid

      CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,            &
                   xs,ys,x2d,y2d,tem1(:,:,1),tem2d_eta,               &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,  &
                   tem4)

      ! 2. Pack and write GRIB message

      CALL mkipds_212(02,102,00,0,p2time,0,0,ipds)
      CALL wrtgb(nx_eta,ny_eta,100*tem2d_eta,nchw,wdlen,grbpkbit,         &
                 ipds,igds, mtem,mitem1,mitem2)
!
!
!-----------------------------------------------------------------------
!
!  Calculate stability indices.
!  Use level k=2 as the "surface".
!
!-----------------------------------------------------------------------
!
!      imid=(nx/2)+1
!      jmid=(ny/2)+1
      CALL arps_be(nx,ny,nz,                                            &
           pprt,zps_km,tem2,qvprt,                                      &
           lcl,lfc,el,twdf,li,pbe,mbe,nbe,tcap,                         &
           wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,                               &
           wrk7,wrk8,wrk9,wrk10,wrk11,wrk12,tem2d1)

!      PRINT *, ' Sample stability values: '
!      PRINT *, '   lcl,lfc : ',lcl(imid,jmid),lfc(imid,jmid)
!      PRINT *, '   el, twdf: ',el(imid,jmid),twdf(imid,jmid)
!      PRINT *, '   li, pbe : ',li(imid,jmid),pbe(imid,jmid)
!      PRINT *, '   mbe, nbe: ',mbe(imid,jmid),nbe(imid,jmid)
!      PRINT *, '   tcap    : ',tcap(imid,jmid)
!
      CALL calcshr(nx,ny,nz,x,y,zp,mf2d,                                &
          pprt,tem2,uprt,vprt,pbe,                                      &
          shr37,ustrm,vstrm,srlfl,srmfl,heli,brn,brnu,blcon,            &
          tem2d1,tem2d2,tem2d3,tem2d4,tem4)

!      PRINT *, ' Sample shear values: '
!      PRINT *, ' shr37,ustrm,vstrm: ',                                  &
!          shr37(imid,jmid),ustrm(imid,jmid),vstrm(imid,jmid)
!      PRINT *, ' srlfl,srmfl,heli: ',                                   &
!          srlfl(imid,jmid),srmfl(imid,jmid),heli(imid,jmid)
!      PRINT *, ' brn,brnu,blcon: ',                                     &
!          brn(imid,jmid),brnu(imid,jmid),blcon(imid,jmid)


!  process near-sfc LI
       PRINT *, ' Storing near-sfc LI'       
       CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,             &
                   xs,ys,x2d,y2d,li,tem2d_eta,                          &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)

       CALL mkipds_212(24,01,00,0,p2time,0,0,ipds)
       CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,          &
                  ipds,igds, mtem,mitem1,mitem2)
 
!  process near-sfc CAPE
       PRINT *, ' Storing near-sfc CAPE'
       CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,             &
                   xs,ys,x2d,y2d,pbe,tem2d_eta,                         &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)

       CALL mkipds_212(157,01,00,0,p2time,0,0,ipds)
       CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,          &
                  ipds,igds, mtem,mitem1,mitem2)

!  process near-sfc CIN      
       PRINT *, ' Storing near-sfc CIN'
 
       CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,             &
                   xs,ys,x2d,y2d,nbe,tem2d_eta,                         &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)

       CALL mkipds_212(156,01,00,0,p2time,0,0,ipds)
       CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,          &
                   ipds,igds, mtem,mitem1,mitem2)

!  process near-sfc Helicity
       PRINT *, ' Storing helicity'
       CALL mkarps2d(nx,ny,nx_eta,ny_eta, iorder,iloc,jloc,             &
                   xs,ys,x2d,y2d,heli,tem2d_eta,                        &
                   dxfld, dyfld,rdxfld, rdyfld,slopey,alphay, betay,    &
                   tem4)

       CALL mkipds_212(190,01,00,0,p2time,0,0,ipds)
       CALL wrtgb(nx_eta,ny_eta,tem2d_eta,nchw,wdlen,grbpkbit,          &
                  ipds,igds, mtem,mitem1,mitem2)

!-----------------------------------------------------------------------
       CLOSE(UNIT=nchw)
       CALL retunit(nchw)      
!
!-----------------------------------------------------------------------
!
!  Create ready file, indicating history dump writing is complete
!
!-----------------------------------------------------------------------
!
       IF( readyfl == 1 ) THEN
         WRITE (filnamr,'(a)') trim(filname) // "_ready"
         CALL getunit( nchoutr )
         OPEN (UNIT=nchoutr,FILE=trim(filnamr))
         WRITE (nchoutr,'(a)') trim(filname)
         CLOSE (nchoutr)
         CALL retunit (nchoutr )
       END IF

    END IF

  END DO


  STOP
END PROGRAM arps2eta212

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

SUBROUTINE extrph(nx,ny,nz,zps,t,pr,iprlvl,hgtcdf) 3
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Extrapolate height by using a std atmos lapse rate
!  below the last physical level.  Above the domain,
!  assume a constant temperature above 300 mb, otherwise
!  use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: zps(nx,ny,nz)
  REAL :: t(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: hgtcdf(nx,ny)
!
  INCLUDE 'phycst.inc'
!
  REAL :: gamma,rddg,const
  PARAMETER ( gamma = 0.0065,    & ! degrees/m  lapse rate
          rddg  = (rd/g),                                               &
              const = (rd*gamma/g) )
!
  INTEGER :: i,j
  REAL :: prcdf
!
  prcdf=100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf < pr(i,j,nz-1)) THEN
        IF(pr(i,j,nz-1) <= 30000.) THEN
          hgtcdf(i,j)=zps(i,j,nz-1) +                                   &
              rddg*t(i,j,nz-1)*ALOG(pr(i,j,nz-1)/prcdf)
        ELSE
          hgtcdf(i,j)=zps(i,j,nz-1) + (t(i,j,nz-1)/gamma)*              &
                      (1.-(prcdf/pr(i,j,nz-1))**const)
        END IF
      ELSE IF(prcdf >= pr(i,j,2)) THEN
        hgtcdf(i,j)=zps(i,j,2) + (t(i,j,2)/gamma)*                      &
               (1.-(prcdf/pr(i,j,2))**const)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrph
!


SUBROUTINE extrpt(nx,ny,nz,t,pr,iprlvl,tcdf) 5
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE EXTRPT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Extrapolate temperature by using a std atmos lapse rate
!  below the last physical level.  Above the domain,
!  assume a constant temperature above 300 mb, otherwise
!  use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  REAL :: t(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: tcdf(nx,ny)

  INCLUDE 'phycst.inc'

  REAL, PARAMETER :: gamma = 0.0065         ! degrees/m  lapse rate
  REAL, PARAMETER :: const = (rd*gamma/g)

  INTEGER :: i,j
  REAL :: prcdf

  prcdf = 100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf <= pr(i,j,nz-1)) THEN
        IF(pr(i,j,nz-1) <= 30000.) THEN
          tcdf(i,j)=t(i,j,nz-1)
        ELSE
          tcdf(i,j)=t(i,j,nz-1)*                                        &
                    ((prcdf/pr(i,j,nz-1))**const)
        END IF
      ELSE IF(prcdf >= pr(i,j,2)) THEN
        tcdf(i,j)=t(i,j,2)*                                             &
                    ((prcdf/pr(i,j,2))**const)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrpt
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE EXTRPQ                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE extrpq(nx,ny,nz,q,pr,iprlvl,qcdf) 17
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign a zero-gradient value to scalars below ground.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: q(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: qcdf(nx,ny)
!
  INTEGER :: i,j
  REAL :: prcdf
!
  prcdf=100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf < pr(i,j,nz-1)) THEN
        qcdf(i,j)=q(i,j,nz-1)
      ELSE IF(prcdf > pr(i,j,2)) THEN
        qcdf(i,j)=q(i,j,2)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrpq

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


SUBROUTINE extrpuv(nx,ny,nz,us,vs,pr,iprlvl,ucdf,vcdf) 3

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign a constant value to sub-terrainian winds
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: us(nx,ny,nz)
  REAL :: vs(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: ucdf(nx,ny)
  REAL :: vcdf(nx,ny)
!
  INTEGER :: i,j
  REAL :: prcdf
!
  prcdf=100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf < pr(i,j,nz-1)) THEN
        ucdf(i,j)=us(i,j,nz-1)
        vcdf(i,j)=vs(i,j,nz-1)
      ELSE IF(prcdf > pr(i,j,2)) THEN
        ucdf(i,j)=us(i,j,2)
        vcdf(i,j)=vs(i,j,2)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrpuv
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINES cal_avor                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Calculate absolute Vort*10^5 (1/s)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Min Zou
!  3/2/98
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------


SUBROUTINE cal_avor(tem9,u,v,x,y,nx,ny,nz,mode,flagsin,omega,           & 2,4
           sinlat,tem1,tem2,tem3)

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: tem9(nx,ny,nz),sinlat(nx,ny)
  REAL :: u(nx,ny,nz), v(nx,ny,nz)
  REAL :: x(nx), y(ny)
  REAL :: tem1(nx), tem2(ny), tem3(nx,ny)
  REAL :: omega
  INTEGER :: mode,flagsin
  
  INTEGER :: i,j,k
  REAL :: tmp1
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=2,nz-2
    DO j=2,ny-2
      DO i=2,nx-2
        tem9(i,j,k)= 1.0E5*(                                            &
               (v(i+1,j,k)-v(i-1,j,k)+v(i+1,j+1,k)-v(i-1,j+1,k))/       &
               (4*(x(i+1)-x(i)))-                                       &
               (u(i,j+1,k)-u(i,j-1,k)+u(i+1,j+1,k)-u(i+1,j-1,k))/       &
               (4*(y(j+1)-y(j))) )
      END DO
    END DO
  END DO

  DO j=2,ny-2
    DO i=2,nx-2
      tem9(i,j,   1)=tem9(i,j,   2)
      tem9(i,j,nz-1)=tem9(i,j,nz-2)
    END DO
  END DO

  DO k=1,nz-1
    DO j=2,ny-2
      tem9(   1,j,k)=tem9(   2,j,k)
      tem9(nx-1,j,k)=tem9(nx-2,j,k)
    END DO
  END DO

  DO k=1,nz-1
    DO i=1,nx-1
      tem9(i,   1,k)=tem9(i,   2,k)
      tem9(i,ny-1,k)=tem9(i,ny-2,k)
    END DO
  END DO

  IF(mode == 1) THEN
    IF( flagsin == 0) THEN
      CALL gtsinlat_ext(nx,ny,x,y,sinlat,tem1,tem2, tem3)
      tmp1 = 2.0* omega
      DO j=1,ny
        DO i=1,nx
          sinlat(i,j) = tmp1 * sinlat(i,j)*1.0E5
        END DO
      END DO
!      flagsin=1
    END IF
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem9(i,j,k) = tem9(i,j,k) + sinlat(i,j)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE cal_avor

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


SUBROUTINE gtsinlat_ext(nx,ny, x,y, sinlat, xs, ys, temxy) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the sin of the lattitude of each grid point, to be used
!  in the calculation of latitude-dependent Coriolis parameters.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Ming Xue
!  3/21/95.
!
!  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)
!    x        x-coordinate of grid points in computational space (m)
!    y        y-coordinate of grid points in computational space (m)
!
!  OUTPUT:
!
!    sinlat   Sin of latitude at each grid point
!
!  WORK ARRAYS:
!
!    xs       x-coordinate at scalar points.
!    ys       y-coordinate at scalar points.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny          ! Dimensions of model grid.
  REAL :: x     (nx)        ! The x-coord. of the u points.
  REAL :: y     (ny)        ! The y-coord. of the v points.

  REAL :: sinlat(nx,ny)     ! Sin of latitude at each grid point

  REAL :: xs    (nx)        ! x-coord. of scalar points.
  REAL :: ys    (ny)        ! y-coord. of scalar points.

  REAL :: temxy (nx,ny)     ! Work array.

  REAL :: r2d
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  DO i=1,nx-1
!    xs(i) = (x(i)+x(i+1))*0.5
!  END DO
!  xs(nx) = -xs(nx-2)+2.0*xs(nx-1)
!
!  DO j=1,ny-1
!    ys(j) = (y(j)+y(j+1))*0.5
!  END DO
!  ys(ny) = -ys(ny-2)+2.0*ys(ny-1)
!
  CALL xytoll(nx,ny,xs,ys,sinlat,temxy)

  r2d = ATAN(1.0)*4.0/180.0

  DO j=1,ny
    DO i=1,nx
      sinlat(i,j) = SIN( r2d * sinlat(i,j) )
    END DO
  END DO

  RETURN
END SUBROUTINE gtsinlat_ext

!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%                                                              %%
!%%  Subroutines for grib dumping                                %%
!%%                                                              %%
!%%  The follwing subroutines were all rewritten from those      %%
!%%  in src/arps/gribio3d.f90 specifically for ETA #212 grid     %%
!%%                                                              %%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE MKIPDS                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE mkipds_212(parno,lvltyp,vlvl,timer,p1time,p2time,scaling, ipds) 28
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Make integer product data array ipds
!
!-----------------------------------------------------------------------
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    parno    Indicator of parameter (see Code table 2)
!    lvltyp   Indicator of type of level (see Code table 3)
!    vlvl     Height, pressure, etc. of levels (see Code table 3)
!    timer    Time range indicator (eithor 10 or 4)
!    p1time
!    p2time   Period of time (number of time units)
!             ARPS forcast time in minutes 
!    scaling  Units decimal scale factor 
!
!  OUTPUT:
!
!    ipds     IPDS array
!
!  WORK ARRAY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: parno,lvltyp,vlvl,timer,p1time,p2time,scaling

  INTEGER, INTENT(OUT) :: ipds(25)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ipds(1) = 28          ! number of bytes in PDS section
  ipds(2) = 2           ! parameter table version number, unknown
            !%% ipds(19) = 2   !? (19)  - VERSION NR OF PARAMETER TABLE

  ipds(3) = 0           ! ID number of originating center. CAPS?
  ipds(4) = 0           ! ID number of model. ARPS/CAPS?
  ipds(5) = 212         ! grid ID, 255 for self-defined GDS
  ipds(6) = 1           ! GDS flag, 1 = GDS section included
            !%%  ipds(1) =  07  ! (1)   - ID OF CENTER
            !%%  ipds(2) = 255  ! (2)   - GENERATING PROCESS ID NUMBER
            !%%  ipds(3) = 212  ! (3)   - GRID DEFINITION
            !%%  ipds(4) = 128  ! (4)   - GDS/BMS FLAG (RIGHT ADJ COPY OF OCTET 8)
                                !         128   GDS included
                                !          64   BMS included
                                !         192   Both included 
  ipds(7) = 0           ! BMS flag, 0 = no BMS section included

!%%%
!%%% After called this subroutine, ipds(8), ipds(9) and ipds(11)
!%%% need to be set for each variables specifically.
!%%%
  ipds(8) = parno       ! indicator of parameter and unit
                        ! (table 2), depends on variables
  ipds(9) = lvltyp      ! indicator of type of level (table 3)
                        ! depends on which variable
                        ! 100 for isobar surface
                        ! 103 for z coordinates in meters
                        ! 111 for soil layers in centimeters
  ipds(10) = 0          ! value 1 of level, N/A
  ipds(11) = vlvl       ! value 2 of level (value of level)
  ! both use two octets: ipds(10) & ipds(11)
           !%% ipds(5) = 011   ! (5)   - INDICATOR OF PARAMETER
           !%% ipds(6) = 100   ! (6)   - TYPE OF LEVEL
           !%% ipds(7) = 1000  ! (7)   - HEIGHT/PRESSURE , ETC 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(8) = year-2000  ! (8)   - YEAR INCLUDING (CENTURY-1)
            !%% ipds(9) = month      ! (9)   - MONTH OF YEAR
            !%% ipds(10)= day        ! (10)  - DAY OF MONTH
            !%% ipds(11)= hour       ! (11)  - HOUR OF DAY
            !%% ipds(12)= minute     ! (12)  - MINUTE OF HOUR
            !%% ipds(21)= 21         ! (21)  - CENTURY OF REFERENCE TIME OF DATA

  ipds(17) = 1          ! forecast time unit, minutes
            !%% ipds(13) = 254  !  (13)  - INDICATOR OF FORECAST TIME UNIT
                                !       254 seconds
                                !         0 minutes
                                !         1 hours
  ipds(18) = p1time     ! forecast time P1, or current time, curtim
  ipds(19) = p2time     ! forecast time P2, N/A
  ipds(20) = timer      ! time range indicator,
                        ! 10 = use two octets from ipds(18)
                        !  0 = Analysis at initial time
                        !  4 = accumulated from p1 to p2
            !%% ipds(14) = 0    !  (14)  - TIME RANGE 1
            !%% ipds(15) = 0    !  (15)  - TIME RANGE 2
            !%% ipds(16) = 1    !  (16)  - TIME RANGE FLAG
  ipds(21) = 0          ! number include in average, N/A
  ipds(22) = 0          ! number missing from average, N/A
            !%% ipds(17) = 0    !  (17)  - NUMBER INCLUDED IN AVERAGE
            !%% ipds(20) = 0    !  (20)  - NR MISSING FROM AVERAGE/ACCUMULATION
  ipds(24) = 0          ! subcenter ID number
  ipds(25) = scaling    ! scaling power of 10,
                        ! depends on which variable
            !%% ipds(23) = 255  !  (23)  - SUBCENTER NUMBER
            !%% ipds(22) = 2    !? (22)  - UNITS DECIMAL SCALE FACTOR

  

! parameters below are those appeared in W3LIB

  !%% ipds(18) = 1   !? (18)  - VERSION NR OF GRIB SPECIFICATION
  !%% ipds(24) = 0   !  (24)  - PDS BYTE 29, FOR NMC ENSEMBLE PRODUCTS
                     !             128 IF FORECAST FIELD ERROR
                     !             64 IF BIAS CORRECTED FCST FIELD
                     !             32 IF SMOOTHED FIELD
                     !      WARNING: CAN BE COMBINATION OF MORE THAN 1
  !%% ipds(25) = 0   !  (25)  - PDS BYTE 30, NOT USED

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


SUBROUTINE mkigds_212(nx,ny,nz,igds) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Make integer grid description array IGDS.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: 
!
!  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:
!
!    igds     IPDS array
!
!  WORK ARRAY:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: nx,ny,nz ! Number of grid points in 3 directions

  INTEGER, INTENT(OUT) :: igds(25)

!-----------------------------------------------------------------------
!
!  ETA #212 Grid parameters
!
!-----------------------------------------------------------------------
  REAL, PARAMETER :: swlat = 12.19, swlon = -133.459
  REAL, PARAMETER :: trulat = 25.0, trulon = -95.0
  REAL, PARAMETER :: dx = 40635.25, dy = 40635.25  ! Meters
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

    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) = NINT(dx)
    igds(11) = NINT(dy)

    igds(12) = 0        ! northern pole on plane

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

    igds(19:25) = 0       ! unused

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


SUBROUTINE wrtgb(nx,ny,fvar,nchanl,wdlen,grbpkbit,              & 28,1
           ipds,igds,tem1,item1,item2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write a 2D float array, fvar, into file pointer nchanl 
!  which points to a GRIB file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!
!  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)
!
!    fvar     Floating array to be written into GRIB file
!
!    nchanl   FILE pointer of GRIB stream file which was opened by
!             a C program, GOPEN
!
!    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).
!    ibdshd   BDS header
!
!  WORK ARRAY:
!
!    pds      PDS (GRIB Section 1)
!    gds      GDS (GRIB Section 3)
!
!    bds      BDS (GRIB Section 4)
!    ibds     Identical to BDS
!
!    nbufsz   Size of buffer of a GRIB message array
!    mgrib    Working buffer of GRIB message array
!
!
!    tem1     Temporary work array.
!    item1    Integer temporary work array.
!    item2    Integer temporary work array.
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  REAL :: fvar(nx,ny)          ! 2D Floating array

  INTEGER :: nchanl            ! FILE pointer indicates the GRIB file
  INTEGER :: wdlen             ! Length of machine word
  INTEGER :: grbpkbit          ! Number of bits in packing GRIB data

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

  REAL :: tem1 (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

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

  INTEGER, PARAMETER :: nbufsz = 800000     ! Size of GRIB message array
  INTEGER :: ibds(nbufsz/4)          ! Identical to BDS
  CHARACTER (LEN=1) :: bds(nbufsz)   ! BDS
  EQUIVALENCE (bds, ibds)

  CHARACTER (LEN=1) :: mgrib(nbufsz) ! Working buffer

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
  INTEGER :: itype             ! Data type: 0 - floating, 1 - integer
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  npts = nx*ny

  itype = 0
  item1 = 0
  ibdshd = 0

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

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

  RETURN
END SUBROUTINE wrtgb
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  FUNCTION compute_density            ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

FUNCTION compute_density(t_k, p_pa) RESULT(rho)

! Computes density (kg/m{-3}) given temperature (K) and pressure (Pa)
  IMPLICIT NONE
  REAL, INTENT(IN) :: t_k
  REAL, INTENT(IN) :: p_pa
  REAL             :: rho

  REAL, PARAMETER  :: Rd =  287.04  ! J deg{-1} kg{-1}
                                    ! Gas constant for dry air

  rho = p_pa / (Rd * t_k)

  RETURN
END FUNCTION compute_density