PROGRAM arps2ncdf,106
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM ARSP2NCDF                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Converts ARPS history files to a netCDF 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.  Sets up variables needed for LDADS/AWIPS D2d display
!  of data.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  02/19/1999
!
!  MODIFICATION HISTORY:
!  08/13/1999
!  Fixed MSLP calculation (J. Levit, K. Brewster)
!
!  1 June 2002 Eric Kemp
!  Soil variable updates.
!
!  11/22/2003 (K. Brewster)
!  Added logic to allow appending new time levels to existing file.
!  Includes additions to input file.  Corrected itim1970 constant.
!
!  12/02/2003 (T. Oram, NWS SFMG)
!  Modifications for AWIPS compatibility plus comments to help
!  myself understand what is going on
!  a.   This program works with the arps.cdl netCDF cdl file that is used
!       for the AWIPS localization process.  The arps.cdl file is new.
!  b.   Changed the tracking of the current and maximum number of files that
!       can be appended to a netCDF file using the record and n_valtimes
!       dimensions in the netCDF file.  
!  c.   modified arps2ncdf.input for new variables and modified write of 
!       1-d variables (valtime, reftime, valtimeMINUSreftime, etc) to only
!       be written once for netcdf append files.  This assumes that all forecast
!       times are known when the first history file is read and converted to
!       netCDF.  Plan to move the number of pressure levels and the pressure level
!       array into arps2ncdf.input for run-time determination rather than
!       netcdf.inc for compile-time determination.
!  d.   The left most dimension of the inventory variables also was 
!       changed to n_valtimes rather than record.  AWIPS chokes 
!       without a full inventory array for all valid times.  The new dimension of
!       the inventory array matches the AWIPS convention.
!  e.   added surface level to the 3-dimensional array to match the awips convention
!  f.   added cw, cice, and tke to 3-dimensional array (tdo, 02/02/2004)
!  g.   removed 2-d fields not included in the dataFieldTable.txt file.
!       also changed reftime to remove the record dimension to match awips (tdo, 02/05/2004)
!
!  02/22/2005 (Y. Wang)
!  Moved the pressure level specification to arps2ncdf.input for run-time
!  determination as planed by T. Oram.
!
!  Cleaned up documents and formats.
!  
!-----------------------------------------------------------------------
!
!  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, PARAMETER  :: nhisfile_max = 200
  INTEGER             :: hinfmt,nhisfile,lengbf,lenfil
  CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'netcdf.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

!
!-----------------------------------------------------------------------
!
!  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
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)
  REAL, ALLOCATABLE :: tem4(:,:,:)
  REAL, ALLOCATABLE :: tem2d1(:,:)
  REAL, ALLOCATABLE :: tem2d2(:,:)
  REAL, ALLOCATABLE :: tem2d3(:,:)
  REAL, ALLOCATABLE :: tem2d4(:,:)

  REAL, ALLOCATABLE :: wrk1(:),wrk2(:),wrk3(:),wrk4 (:),wrk5 (:),wrk6 (:)
  REAL, ALLOCATABLE :: wrk7(:),wrk8(:),wrk9(:),wrk10(:),wrk11(:),wrk12(:)

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

! J.Case, ENSCO Inc. (9/29/2004)
! Additional variables.
! Note: radar_comp is maximum column reflecitivity and
!       echo_top is maximum height of 18 dBZ echo in a column.

  REAL, ALLOCATABLE :: radar_3d(:,:,:),radar_comp(:,:),echo_top(:,:)

  REAL, ALLOCATABLE :: mf2d(:,:)
  REAL, ALLOCATABLE :: coriolis(:,:)
  REAL, ALLOCATABLE :: spacing(:,:)
!
!-----------------------------------------------------------------------
!
!  Misc ARPS variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nch

  INTEGER :: first18                    ! J.Case, ENSCO Inc. (9/29/2004)

  REAL    :: time
  REAL    :: latnot(2)
  LOGICAL :: init
!
!-----------------------------------------------------------------------
!
!  netCDF dimensions
!
!  nprlvl       defined in arps2cdf.inc is the number of pressure levels
!  nxcdf        number of points in x direction in netcdf file
!  nycdf        number of points in y direction in netcdf file
!  nrecord      number of records in netcdf file
!  ntotaltime   number of valid times that are to be or can be written into file
!  ncdflvls     number of levels in 3-d grids in file; nprlvl + 1 (surface)
!
!-----------------------------------------------------------------------
!
  INTEGER, PARAMETER :: mxtime = 40
  INTEGER :: nxcdf,nycdf,nrecord,nvaltimes,ntotaltimes
  
  INTEGER, PARAMETER :: nprlvl_max = 128
  INTEGER            :: iprlvl(nprlvl_max)
  INTEGER            :: nprlvl, ncdflvls                     ! = nprlvl+1

  INTEGER, PARAMETER :: n3dvar = 10, n2dvar = 16, nstvar = 3
  INTEGER, PARAMETER :: nvartot = (n3dvar+n2dvar)
                        ! J.Case, ENSCO Inc. (9/29/2004) -- Added 1 new 3D 
                        ! variable (rr) and 2 new 2D variables (cxr, mret)
!
!-----------------------------------------------------------------------
!
!  netCDF variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: istatus,ncid
  INTEGER :: dimidrec,dimidtot,dimidnt,dimidnav
  INTEGER :: dimidnx,dimidny,dimid1,dimidnz,dimchrlvl,dimnamlen
  INTEGER :: vtmrefid,valtimid,reftimid,origid,modelid,ntwrtid

  INTEGER :: vecistr(2)
  INTEGER :: planeistr(4)
  INTEGER :: volistart(4)
  DATA vecistr /1,1/
  DATA planeistr /1,1,1,1/
  DATA volistart /1,1,1,1/

  INTEGER :: planedim(4)
  INTEGER :: voldim(4)
  INTEGER :: lvldim(2)
  INTEGER :: invdim(2)
  INTEGER :: planelen(4)
  INTEGER :: lvllen(2)
  INTEGER :: sfclen(2)
  INTEGER :: vollen(4)
  INTEGER :: invlen(2)
  INTEGER :: ntwrt

  REAL :: vrange(2)
!
!-----------------------------------------------------------------------
!
!  netCDF output variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: v3did(n3dvar)
  INTEGER :: v3dlvlid(n3dvar)
  INTEGER :: v2dinvid(n2dvar)
  INTEGER :: v3dinvid(n3dvar)
  INTEGER :: v2did(n2dvar)
  INTEGER :: vstid(nstvar)
  INTEGER :: v2dlvlid(n2dvar)
  INTEGER :: vtmref(mxtime)
  REAL*8  valtime(mxtime)
  REAL*8  reftime(mxtime)

  CHARACTER (LEN=4) :: v3dnam(n3dvar)
  CHARACTER (LEN=4) :: v2dnam(n2dvar)
  CHARACTER (LEN=14) :: vstnam(nstvar)
  CHARACTER (LEN=10) :: v3dlvl(nprlvl_max)
  CHARACTER (LEN=10) :: v2dlvl
  CHARACTER (LEN=30) :: v3dlngnam(n3dvar)
  CHARACTER (LEN=30) :: v2dlngnam(n2dvar)
  CHARACTER (LEN=30) :: vstlngnam(nstvar)
  CHARACTER (LEN=30) :: v3dunits(n3dvar)
  CHARACTER (LEN=30) :: v2dunits(n2dvar)
  CHARACTER (LEN=30) :: vstunits(nstvar)
  CHARACTER (LEN=nprlvl_max) :: inventory
  CHARACTER (LEN=1) :: inventory_sfc
  CHARACTER (LEN=1) :: inventory_blank
  REAL :: v3dmin(n3dvar)
  REAL :: v2dmin(n2dvar)
  REAL :: v3dmax(n3dvar)
  REAL :: v2dmax(n2dvar)

  CHARACTER (LEN=10) :: v2dlvlnam
  CHARACTER (LEN=10) :: v3dlvlnam
  CHARACTER (LEN=13) :: v2dinv
  CHARACTER (LEN=13) :: v3dinv

  DATA vtmref /mxtime*0/
!
!-----------------------------------------------------------------------
!
!  Variable indexes and descriptors
!
!-----------------------------------------------------------------------
!
  INTEGER :: idxgh,idxrh,idxt,idxuw,idxvw,idxww,idxcw,idxci,idxtk
                                         ! J.Case, ENSCO Inc. (9/29/2004)
  INTEGER :: idxrr
  PARAMETER (idxgh=1,idxrh=2,idxt=3,                                    &
             idxuw=4,idxvw=5,idxww=6,                                   &
             idxcw=7,idxci=8,idxtk=9,idxrr=10)
!
! J.Case, cont. (9/29/2004) -- Added radar reflectivity to 3D variables (rr).

  DATA v3dnam /'gh','rh','t','uw','vw','ww','cw','cice','tke','rr'/
  DATA v3dlngnam /'Geopotential Height',                                &
                  'Relative Humidity',                                  &
                  'Temperature',                                        &
                  'u Wind Component',                                   &
                  'v Wind Component',                                   &
                  'w Wind Component',                                   &
                  'Cloud Water',                                        &
                  'Cloud Ice',                                          &
                  'Turbulent Kinetic Energy',                           &
                  'Radar Reflectivity'/
  DATA v3dunits /'geopotential meters',                                 &
                 'percent',                                             &
                 'degrees K',                                           &
                 'meters/second',                                       &
                 'meters/second',                                       &
                 'meters/second',                                       &
                 'kg/kg',                                               &
                 'kg/kg',                                               &
                 '(meters/second)**2',                                  &
                 'dBZ'/
  DATA v3dmin /    0.,  0.,150.,-300.,-300.,-100.,0.,0.,0.,-200./
  DATA v3dmax /40000.,110.,400., 300., 300., 100.,1.,1.,1000.,200./

  INTEGER :: idxmslp,idxp,idxlgsp,idxcp,idxtp
  INTEGER :: idxdpt
  INTEGER :: idxsh,idxept,idxsli,idxcape,idxcin
  INTEGER :: idxustrm,idxvstrm
  INTEGER :: idxheli

         ! J.Case, ENSCO Inc. (9/29/2004) -- new 2D variables cxr and mret
  INTEGER :: idxcxr,idxmret

  PARAMETER (idxmslp=1,                                                 &
             idxp   =2,                                                 &
             idxlgsp=3,                                                 &
             idxcp  =4,                                                 &
             idxtp  =5,                                                 &
             idxdpt =6,                                                 &
             idxsh  =7,                                                 &
             idxept =8,                                                 &
             idxsli =9,                                                 &
             idxcape=10,                                                &
             idxcin =11,                                                &
             idxustrm=12,                                               &
             idxvstrm=13,                                               &
             idxheli =14,                                               &
             idxcxr  =15,                                               &
             idxmret =16)

!             idxsfu =6,                                                 &
!             idxsfv =7,                                                 &
!             idxsft =8,                                                 &
!             idxlcl =15,                                                &
!             idxlfc =16,                                                &
!             idxel  =17,                                                &
!             idxtwdf =18,                                               &
!             idxtcap =19,                                               &
!             idxshr37=20,                                               &
!             idxsrlfl=23,                                               &
!             idxsrmfl=24,                                               &
!             idxbrn  =26,                                               &
!             idxbrnu =27,                                               &
!             idxblcon=28)


! these 2d variable names must match cdl names in the AWIPS
! dataFieldTable.txt file
!  DATA v2dnam /'mslp','p','lgsp','cp','tp','sfu','sfv',                 &
!               'sft','dpt','sh','ept','sli','cape','cin','lcl',         &
!               'lfc','el','twdf','tcap','sh37','ustm','vstm',           &
!               'srlf','srmf','heli','brn','brnu','blcn'/

!JLC (8/23/2004) -- added cxr,mret to 2D variables.

  DATA v2dnam /'mslp','p','lgsp','cp','tp',                             &
               'dpt','sh','ept','sli','cape','cin',                     &
               'ustm','vstm',                                           &
               'heli','cxr','mret'/


  DATA v2dlngnam /'Mean Sea Level Pressure',                            &
                  'Surface Pressure',                                   &
                  'Grid Scale Precipitation',                           &
                  'Convective Precipitation',                           &
                  'Total Precipitation',                                &
                  'Surface Dew Point',                                  &
                  'Specific Humidity',                                  &
                  'Theta-e',                                            &
                  'Surface Lifted Index',                               &
                  'CAPE',                                               &
                  'CIN',                                                &
                  'Est Storm Motion u comp',                            &
                  'Est Storm Motion v comp',                            &
                  'Storm Rel Helicity',                                 &
                  'Column Max Reflectivity',                            &
                  'Maximum Radar Echo Tops'/

!                  'Surface u Wind Component',                           &
!                  'Surface v Wind Component',                           &
!                  'Surface Temperature',                                &
!                  'Lifting Condensation Lvl',                           &
!                  'Level of Free Convection',                           &
!                  'Equilibrium Level',                                  &
!                  'Max Wet-Bulb Temp Diff',                             &
!                  'Cap Strength',                                       &
!                  '3-7km Shear',                                        &
!                  'Storm Rel Low-Level Flow',                           &
!                  'Storm Rel Mid-Level Flow',                           &
!                  'Bulk Richardson Number',                             &
!                  'BRN Shear u',                                        &
!                  'Bound-Layer Conv x 1000'/

! J.Case, cont. (9/29/2004)

  DATA v2dunits /'Pascals','Pascals','mm','mm','mm',                    &
                 'degrees K','kg/kg',                                   &
                 'degrees K','degrees C','J/kg','J/kg',                 &
                 'meters/second','meters/second',                       &
                 'm2/s2','dBZ','meters'/

  DATA v2dmin / 80000.,50000.,0.,0.,0.,                                 &
                -100.,0.,                                               &
                200.,-100.,0.,0.,                                       &
                -100.,-100.,                                            &
                -900.,0.,0./

  DATA v2dmax / 110000.,120000.,1000.,1000.,1000.,                      &
                100.,1.,                                                &
                500.,100.,10000.,10000.,                                &
                100.,100.,                                              &
                900.,100.,20000./

  INTEGER :: idxtop,idxcor,idxspa
  PARAMETER (idxtop=1,idxcor=2,idxspa=3)
  DATA vstnam /'staticTopo',                                            &
               'staticCoriolis',                                        &
               'staticSpacing'/
  DATA vstlngnam /'Topography',                                         &
                  'Coriolis parameter',                                 &
                  'Grid spacing'/
  DATA vstunits /'meters',                                              &
                 '/second',                                             &
                 'meters'/

  CHARACTER (LEN=20) :: cproj
!
!  The following are the projections supported by LDAD/AWIPS
!  in their ordering/naming convention
!
  CHARACTER (LEN=25) :: chproj(17)
  DATA chproj /'STEREOGRAPHIC',                                         &
               'ORTHOGRAPHIC',                                          &
               'LAMBERT_CONFORMAL',                                     &
               'AZIMUTHAL_EQUAL_AREA',                                  &
               'GNOMONIC',                                              &
               'AZIMUTHAL_EQUIDISTANT',                                 &
               'SATELLITE_VIEW',                                        &
               'CYLINDRICAL_EQUIDISTANT',                               &
               'MERCATOR',                                              &
               'MOLLWEIDE',                                             &
               'NULL',                                                  &
               'NULL',                                                  &
               'NULL',                                                  &
               'NULL',                                                  &
               'NULL',                                                  &
               'AZIMUTH_RANGE',                                         &
               'RADAR_AZ_RAN'/
!
!  kproj is the index in the LDAD/AWIPS map projection list
!  for the mapproj number in the ARPS file
!
  INTEGER :: kproj(3)
  DATA kproj /1,3,9/

  REAL, ALLOCATABLE :: static2d(:,:,:)
  REAL, ALLOCATABLE :: cdf2d(:,:,:,:)
  REAL, ALLOCATABLE :: cdf3d(:,:,:,:)
!
!-----------------------------------------------------------------------
!
!  Constants
!
!-----------------------------------------------------------------------
!
  INTEGER, PARAMETER  :: itim1970=315622800
  CHARACTER(LEN=8),   PARAMETER :: cdldate  = '20031202'
  CHARACTER(LEN=40),  PARAMETER :: depictor = 'ARPS'
  CHARACTER(LEN=132), PARAMETER :: origin   = 'CAPS'
  CHARACTER(LEN=132), PARAMETER :: model    = 'ARPS 5.3'
!
!-----------------------------------------------------------------------
!
!  Namelists
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=132) :: cdfname
  INTEGER :: ipktyp,nbits,ncdapnd,ncdntime,ncdtinterval
  NAMELIST /ncdf_options/ nprlvl,iprlvl,cdfname,ipktyp,nbits,ncdapnd, &
                          ncdntime,ncdtinterval
!
!-----------------------------------------------------------------------
!
!  External functions
!
!-----------------------------------------------------------------------
!
  REAL  :: wmr2td,oe,dpt
  EXTERNAL wmr2td,oe,dpt
!
!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=10) :: nzstr
  REAL    :: latll,lonll,latur,lonur,dxkm,dykm,latdxdy,londxdy,rotate
  INTEGER :: ivar,ls,itime,itimref,iproj,idxtime, ntime

  INTEGER :: i,j,k,klev,ifile,ireturn,imid,jmid,ncmode
  REAL    :: rlnpcdf,r2d
  REAL    :: ctrx,ctry,swx,swy
  REAL    :: gamma,ex1,ex2,rln700,p00,t00
  REAL    :: prmb,tdew
  LOGICAL :: adddefs,ncexist

  REAL, ALLOCATABLE :: p_mb(:,:,:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Initializations
!
!-----------------------------------------------------------------------
!
!
  WRITE(6,'(9(/a))')                                                    &
      '###############################################################',&
      '###############################################################',&
      '####                                                       ####',&
      '####                Welcome to ARPS2NCDF                   ####',&
      '####        A program that reads in history files          ####',&
      '####  generated by ARPS and produces a netCDF format file. ####',&
      '####                                                       ####',&
      '###############################################################',&
      '###############################################################'
!
!-----------------------------------------------------------------------
!
!  Get the names of the input data files.
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  get_input_filenames reads the history file input options 
!  from the namelist then open the file and read some of 
!  the grid dimension information
!
!-----------------------------------------------------------------------
!
  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

!
!-----------------------------------------------------------------------
!
!  Read the netcdf output options from the netcdf namelist options
!       ncdapnd    option to append to netcdf file
!       ncdntime   the number of valid times in the file
!       the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
  ncdapnd=0
  ipktyp=1
  nbits=16
  READ(5,ncdf_options,END=900)

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

!-----------------------------------------------------------------------
! output grid is 3 gridpoints less than history file grid 
! in both x and y directions
!-----------------------------------------------------------------------

  ncdflvls = nprlvl + 1       ! include surface layer

  nxcdf=(nx-3)
  nycdf=(ny-3)

  planelen(1) = nxcdf
  planelen(2) = nycdf
  planelen(3) = 1
  planelen(4) = 1
  lvllen(1)   = 10
  lvllen(2)   = ncdflvls
  sfclen(1)   = 10
  sfclen(2)   = 1
  vollen(1)   = nxcdf
  vollen(2)   = nycdf
  vollen(3)   = ncdflvls
  vollen(4)   = 1
  ncdntime    = max(nhisfile,ncdntime)

!-----------------------------------------------------------------------
!
! now that you've read all the dimensions at execute time
! allocate the memory for the arrays
!
!-----------------------------------------------------------------------
  ALLOCATE(x      (nx))
  ALLOCATE(y      (ny))
  ALLOCATE(z      (nz))
  ALLOCATE(zp     (nx,ny,nz))
  ALLOCATE(zpsoil (nx,ny,nzsoil))

  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(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))

  ALLOCATE(e_mb   (nx,ny,nz))
  ALLOCATE(mix    (nx,ny,nz))
  ALLOCATE(esat_mb(nx,ny,nz))
  ALLOCATE(rh     (nx,ny,nz))
  ALLOCATE(t_dew  (nx,ny,nz))

  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))

  ALLOCATE(tem1(nx,ny,nz))
  ALLOCATE(tem2(nx,ny,nz))
  ALLOCATE(tem3(nx,ny,nz))
  ALLOCATE(tem4(nx,ny,nz))
  ALLOCATE(tem2d1(nx,ny))
  ALLOCATE(tem2d2(nx,ny))
  ALLOCATE(tem2d3(nx,ny))
  ALLOCATE(tem2d4(nx,ny))

  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))

  ALLOCATE(xs(nx))
  ALLOCATE(ys(ny))
  ALLOCATE(zps(nx,ny,nz))
  ALLOCATE(zps_km(nx,ny,nz))

  ! J.Case, ENSCO Inc. (9/29/2004)
  ALLOCATE(radar_3d(nx,ny,nz))
  ALLOCATE(radar_comp(nx,ny))
  ALLOCATE(echo_top(nx,ny))

  ALLOCATE(p_mb(nx,ny,nz))

  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))

  ALLOCATE(mf2d(nx,ny))
  ALLOCATE(coriolis(nx,ny))
  ALLOCATE(spacing(nx,ny))


  ALLOCATE(static2d(nxcdf,nycdf,nstvar))
  ALLOCATE(cdf2d(nxcdf,nycdf,1,n2dvar))
  ALLOCATE(cdf3d(nxcdf,nycdf,ncdflvls,n3dvar))

!-----------------------------------------------------------------------
! initialize a bunch of stuff
!-----------------------------------------------------------------------
  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

  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

  e_mb    =0.0
  mix     =0.0
  esat_mb =0.0
  rh      =0.0
  t_dew   =0.0

  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

  tem1 =0.0
  tem2 =0.0
  tem3 =0.0
  tem4 =0.0
  tem2d1=0.0
  tem2d2=0.0
  tem2d3=0.0
  tem2d4=0.0

  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

  xs=0.0
  ys=0.0
  zps=0.0
  zps_km=0.0

  ! J.Case, ENSCO Inc. (9/29/2004)
  radar_3d=0.0
  radar_comp=0.0
  echo_top=0.0

  p_mb =0.0

  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

  mf2d=0.0
  coriolis=0.0
  spacing=0.0

  init=.false.

  nrecord=0
  nvaltimes=0

  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,nch,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)

!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
    IF( ireturn == 0 ) THEN   ! successful read
!
!===============================================================
!
!  J.Case, ENSCO, Inc. 9/29/2004 -- added subroutine call for 
!                          derived radar reflectivity products.
!
!===============================================================
!

      WRITE(*,*) ' '
      WRITE (*,*) 'Deriving radar reflectivity.'
      CALL REFLEC (nx,ny,nz, rhobar, qr, qs, qh, radar_3d)

      DO j=1,ny-1
        DO i=1,nx-1
          radar_comp(i,j) = 0.0
          echo_top(i,j)   = 0.0
          first18         = 1

          DO k=1,nz-1
            IF (radar_3d(i,j,k) > radar_comp(i,j)) THEN
              radar_comp(i,j) = radar_3d(i,j,k)
            ENDIF
          END DO

          DO k=nz-1,1,-1
            IF (radar_3d(i,j,k) >= 18.0 .AND. first18 == 1) THEN
              echo_top(i,j) = zp(i,j,k)
              first18       = 0
            ENDIF
          END DO

        END DO
      END DO

      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
!
!-----------------------------------------------------------------------
!
!  Global Attributes
!
!-----------------------------------------------------------------------
!
        latnot(1)=trulat1
        latnot(2)=trulat2
        CALL setmapr( mapproj, 1.0 , latnot , trulon)

        CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
        dx=x(3)-x(2)
        dy=y(3)-y(2)
        dxkm=0.001*dx
        dykm=0.001*dy
        swx = ctrx - (FLOAT(nx-3)/2.) * dx
        swy = ctry - (FLOAT(ny-3)/2.) * dy
        CALL setorig( 1, swx, swy)

        CALL xytoll(1,1,xs(2),ys(2),latll,lonll)
        CALL xytoll(1,1,xs(nx-2),ys(ny-2),latur,lonur)

        IF(mapproj > 0) THEN
          iproj=kproj(mapproj)
          cproj=chproj(iproj)
        ELSE
          iproj=0
          cproj='NONE'
        END IF

        IF(mapproj == 2) THEN
          latdxdy=0.5*(trulat1+trulat2)
        ELSE
          latdxdy=trulat1
        END IF
        londxdy=ctrlon

        CALL xytoll(nx,ny,xs,ys,coriolis,tem1)
        CALL xytomf(nx,ny,xs,ys,mf2d)

        r2d = ACOS(-1.0)/180.0
        DO j=1,ny-1
          DO i=1,nx-1
            coriolis(i,j) = 2.*omega*SIN( r2d * coriolis(i,j) )
            spacing(i,j) = dx*mf2d(i,j)
          END DO
        END DO
!-----------------------------------------------------------------------
!
!  Open netCDF file
!
!-----------------------------------------------------------------------
!
        INQUIRE(FILE=cdfname,EXIST=ncexist)

!-----------------------------------------------------------------------
!       if this is an append and the netCDF files exists
!       then read the record and n_valtimes dimensions from the 
!       netcdf file
!-----------------------------------------------------------------------
        IF(ncdapnd > 0 .AND. ncexist) THEN
          ncmode=NF_WRITE
          adddefs=.false.
          istatus=NF_OPEN(cdfname,ncmode,ncid)
          IF (istatus /= NF_NOERR) CALL handle_err('nf_open',istatus)

          ! i think we want to replace with a read of the 
          ! record dimension to determine the record to write to
          ! and the n_valtimes dimension to determine the maximum
          ! number of forecast times that can be included in the file (tdo)

          istatus=NF_INQ_DIMID(ncid,'record',dimidrec)
          IF (istatus /= NF_NOERR) CALL handle_err('nf_inq_dimid record',istatus)
          istatus=NF_INQ_DIMLEN(ncid,dimidrec,nrecord)

          istatus=NF_INQ_DIMID(ncid,'n_valtimes',dimidnt)
          IF (istatus /= NF_NOERR) CALL handle_err('nf_inq_dimid ntime',istatus)
          istatus=NF_INQ_DIMLEN(ncid,dimidnt,nvaltimes)

          WRITE(6,'(2(a,i5))') 'current record =',nrecord,              &
                ',number of valid times =',nvaltimes,                   &
                ',number of history files =',nhisfile
          ntotaltimes=nrecord+nhisfile
          IF (ntotaltimes > nvaltimes) THEN
            WRITE(6,'(a)') 'Trying to write too many hours into netcdf file'
            STOP
          ENDIF

          ntwrt=nrecord
          ! i think this is the simplest change to preserve Keith Brewsters
          ! append capability while living within the AWIPS netCDF conventions

!          istatus=nf_inq_varid(ncid,'ntwrite',ntwrtid)
!          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus)
!          istatus=nf_get_var_int(ncid,ntwrtid,ntwrt)
!          IF (istatus /= nf_noerr) CALL handle_err('nf_get_var_int',istatus)

        ELSE
          ncmode=NF_CLOBBER
          adddefs=.true.
          istatus=NF_CREATE(cdfname,ncmode,ncid)
          ntwrt=0
          IF (istatus /= NF_NOERR) CALL handle_err('nf_create',istatus)
        END IF
!
!-----------------------------------------------------------------------
!
!  Define dimensions
!
!-----------------------------------------------------------------------
!
        IF(ncdflvls < 10) THEN
          WRITE(nzstr,'(a,i1)') 'levels_',ncdflvls
        ELSE IF(ncdflvls < 100) THEN
          WRITE(nzstr,'(a,i2)') 'levels_',ncdflvls
        ELSE
          WRITE(nzstr,'(a,i3)') 'levels_',ncdflvls
        END IF

        IF( adddefs ) THEN

          istatus=nf_def_dim(ncid,'record',nf_unlimited,dimidrec)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim record',istatus)

          istatus=nf_def_dim(ncid,'n_valtimes',ncdntime,dimidnt)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ncdntime',istatus)

          istatus=nf_def_dim(ncid,'data_variables',nvartot,dimidtot)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nvartot',istatus)

          istatus=nf_def_dim(ncid,'charsPerLevel',10,dimchrlvl)
          IF (istatus /= nf_noerr)                                        &
              CALL handle_err('nf_def_dim charsPerLevel',istatus)

          istatus=nf_def_dim(ncid,'namelen',132,dimnamlen)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim namelen',istatus)

          istatus=nf_def_dim(ncid,'x',nxcdf,dimidnx)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nx',istatus)

          istatus=nf_def_dim(ncid,'y',nycdf,dimidny)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ny',istatus)

          istatus=nf_def_dim(ncid,'levels_1',1,dimid1)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim 1',istatus)

          istatus=nf_def_dim(ncid,nzstr,ncdflvls,dimidnz)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim levels',istatus)

          istatus=nf_def_dim(ncid,'nav',1,dimidnav)
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim nav',istatus)
!
!       The following arrays are used to create the dimensions of the variables
!       that are defined next.  The arrays consist of dimension ids (from above)
!       When listed by the ncdump utility, the 1st dimension (i.e. arraydim(1)) is
!       the right most dimension listed while the last dimension is the left most
!       dimension.  The order of dimensions here is set to match the AWIPS convention
!
        planedim(1)=dimidnx
        planedim(2)=dimidny
        planedim(3)=dimid1
        planedim(4)=dimidrec
        voldim(1)=dimidnx
        voldim(2)=dimidny
        voldim(3)=dimidnz
        voldim(4)=dimidrec
        lvldim(1)=dimchrlvl
        lvldim(2)=dimidnz
        invdim(1)=dimidnz
        invdim(2)=dimidnt
!
!-----------------------------------------------------------------------
!
!  Define 3-d variables and attributes
!
!-----------------------------------------------------------------------
!
        DO ivar=1,n3dvar
          vrange(1)=v3dmin(ivar)
          vrange(2)=v3dmax(ivar)

          istatus=nf_def_var(ncid,v3dnam(ivar),nf_float,4,              &
                             voldim,v3did(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var 3d',istatus)

          ls=LEN(v3dlngnam(ivar))
          CALL strlnth(v3dlngnam(ivar),ls)
          istatus=nf_put_att_text(ncid,v3did(ivar),'long_name',         &
                          ls,v3dlngnam(ivar)(1:ls))
          ls=LEN(v3dunits(ivar))
          CALL strlnth(v3dunits(ivar),ls)
          istatus=nf_put_att_text(ncid,v3did(ivar),'units',             &
                          ls,v3dunits(ivar)(1:ls))
          istatus=nf_put_att_real(ncid,v3did(ivar),'valid_range',       &
                          nf_float,2,vrange)
          istatus=nf_put_att_real(ncid,v3did(ivar),'_FillValue',        &
                          nf_float,1,-99999.)
          istatus=nf_put_att_int(ncid,v3did(ivar),'_n3D',               &
                          nf_int,1,ncdflvls)
          istatus=nf_put_att_text(ncid,v3did(ivar),'levels',            &
                          10,'MB and SFC')

          ls=LEN(v3dnam(ivar))
          CALL strlnth(v3dnam(ivar),ls)
          WRITE(v3dlvlnam,'(a,a)') v3dnam(ivar)(1:ls),'Levels'
          istatus=nf_def_var(ncid,v3dlvlnam,nf_char,2,                  &
                             lvldim,v3dlvlid(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var lvl 3d',istatus)
!
!         create the inventory variable name by appending Inventory to variable
!         then define the inventory variable
!
          WRITE(v3dinv,'(a,a)') v3dnam(ivar)(1:ls),'Inventory'
          istatus=nf_def_var(ncid,v3dinv,nf_char,2,                     &
                             invdim,v3dinvid(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var inv 3d',istatus)

        END DO
!
!-----------------------------------------------------------------------
!
!  Define 2-d variables and attributes
!
!-----------------------------------------------------------------------
!
        invdim(1)=dimid1
        invdim(2)=dimidnt
        lvldim(2)=dimid1
        DO ivar=1,n2dvar
          vrange(1)=v2dmin(ivar)
          vrange(2)=v2dmax(ivar)
          istatus=nf_def_var(ncid,v2dnam(ivar),nf_float,4,              &
                             planedim,v2did(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var 2d',istatus)
          ls=LEN(v2dlngnam(ivar))
          CALL strlnth(v2dlngnam(ivar),ls)
          istatus=nf_put_att_text(ncid,v2did(ivar),'long_name',         &
                                  ls,v2dlngnam(ivar)(1:ls))
          ls=LEN(v2dunits(ivar))
          CALL strlnth(v2dunits(ivar),ls)
          istatus=nf_put_att_text(ncid,v2did(ivar),'units',             &
                                  ls,v2dunits(ivar)(1:ls))
          istatus=nf_put_att_real(ncid,v2did(ivar),'valid_range',       &
                                  nf_float,2,vrange)
          istatus=nf_put_att_real(ncid,v2did(ivar),'_FillValue',        &
                          nf_float,1,-99999.)
          istatus=nf_put_att_int(ncid,v2did(ivar),'_n3D',               &
                          nf_int,1,0)
          IF(ivar == idxmslp) THEN
            istatus=nf_put_att_text(ncid,v2did(ivar),'levels',          &
                          3,'MSL')
          ELSE
            istatus=nf_put_att_text(ncid,v2did(ivar),'levels',          &
                          3,'SFC')
          END IF

          ls=LEN(v2dnam(ivar))
          CALL strlnth(v2dnam(ivar),ls)
          WRITE(v2dlvlnam,'(a,a)') v2dnam(ivar)(1:ls),'Levels'
          istatus=nf_def_var(ncid,v2dlvlnam,nf_char,2,                  &
                             lvldim,v2dlvlid(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var lvl 2d',istatus)
!
!         create the inventory variable name by appending Inventory to variable
!         then define the inventory variable
!
          WRITE(v2dinv,'(a,a)') v2dnam(ivar)(1:ls),'Inventory'
          istatus=nf_def_var(ncid,v2dinv,nf_char,2,                     &
                             invdim,v2dinvid(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var inv 2d',istatus)
        END DO
!
!-----------------------------------------------------------------------
!
!  Define 1-d variables
!
!-----------------------------------------------------------------------
!
        istatus=nf_def_var(ncid,'valtimeMINUSreftime',nf_int,1,         &
                           dimidnt,vtmrefid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var vtmref',istatus)
        istatus=nf_put_att_text(ncid,vtmrefid,'units',                  &
                                7,'seconds')

        istatus=nf_def_var(ncid,'reftime',nf_double,0,                  &
                           0,reftimid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var reftim',istatus)
        istatus=nf_put_att_text(ncid,reftimid,'units',                  &
                   33,'seconds since 1970-1-1 00:00:00.0')

        istatus=nf_def_var(ncid,'valtime',nf_double,1,                  &
                           dimidnt,valtimid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var valtime',istatus)
        istatus=nf_put_att_text(ncid,vtmrefid,'units',                  &
                   33,'seconds since 1970-1-1 00:00:00.0')

        istatus=nf_def_var(ncid,'origin',nf_char,1,                     &
                           dimnamlen,origid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var orig',istatus)

        istatus=nf_def_var(ncid,'model',nf_char,1,                      &
                           dimnamlen,modelid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var model',istatus)

!       commented out the following definition of ntwrt variable
!       probably should be a global attribute (tdo)
!        istatus=nf_def_var(ncid,'ntwrite',nf_int,0,0,ntwrtid)
!        IF (istatus /= nf_noerr) CALL handle_err('nf_def_var ntwrite',istatus)

!
!  Define static variables
!
        DO ivar=1,nstvar
          istatus=nf_def_var(ncid,vstnam(ivar),nf_float,2,              &
                             planedim,vstid(ivar))
          IF (istatus /= nf_noerr) CALL handle_err('nf_def_var static',istatus)
          ls=LEN(vstlngnam(ivar))
          CALL strlnth(vstlngnam(ivar),ls)
          istatus=nf_put_att_text(ncid,vstid(ivar),'long_name',         &
                                  ls,vstlngnam(ivar)(1:ls))
          ls=LEN(vstunits(ivar))
          CALL strlnth(vstunits(ivar),ls)
          istatus=nf_put_att_text(ncid,vstid(ivar),'units',             &
                                  ls,vstunits(ivar)(1:ls))
        END DO

        WRITE(6,'(a,a)') '  cdl date ',cdldate

        istatus=nf_put_att_text(ncid,nf_global,'cdlDate',               &
                                8,cdldate)
        ls=LEN(depictor)
        CALL strlnth(depictor,ls)
        istatus=nf_put_att_text(ncid,nf_global,'depictorName',          &
                                ls,depictor)
        istatus=nf_put_att_int(ncid,nf_global,'projIndex',              &
                                nf_int,1,iproj)
        ls=LEN(cproj)
        CALL strlnth(cproj,ls)
        istatus=nf_put_att_text(ncid,nf_global,'projName',              &
                                ls,cproj)
        istatus=nf_put_att_real(ncid,nf_global,'centralLat',            &
                                nf_float,1,ctrlat)
        istatus=nf_put_att_real(ncid,nf_global,'centralLon',            &
                                nf_float,1,ctrlon)
        rotate=ctrlon-trulon
        istatus=nf_put_att_real(ncid,nf_global,'rotation',              &
                                nf_float,1,rotate)
        istatus=nf_put_att_real(ncid,nf_global,'lat00',                 &
                                nf_float,1,latll)
        istatus=nf_put_att_real(ncid,nf_global,'lon00',                 &
                                nf_float,1,lonll)
        istatus=nf_put_att_real(ncid,nf_global,'latNxNy',               &
                                nf_float,1,latur)
        istatus=nf_put_att_real(ncid,nf_global,'lonNxNy',               &
                                nf_float,1,lonur)
        istatus=nf_put_att_real(ncid,nf_global,'dxKm',                  &
                                nf_float,1,dxkm)
        istatus=nf_put_att_real(ncid,nf_global,'dyKm',                  &
                                nf_float,1,dykm)
        istatus=nf_put_att_real(ncid,nf_global,'latDxDy',               &
                                nf_float,1,latdxdy)
        istatus=nf_put_att_real(ncid,nf_global,'lonDxDy',               &
                                nf_float,1,londxdy)

!
!-----------------------------------------------------------------------
!
!  End define mode
!
!-----------------------------------------------------------------------
!
          istatus=nf_enddef(ncid)
          IF (istatus /= nf_noerr) CALL handle_err('nf_enddef',istatus)

        ELSE  ! get variable ids from existing dataset

!
!  Get dimension ids
!  The inquery for record and n_valtimes is redundant now, but ok so won't touch 
!
          istatus=nf_inq_dimid(ncid,'record',dimidrec)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid record',istatus)

          istatus=nf_inq_dimid(ncid,'n_valtimes',dimidnt)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid ntime',istatus)

          istatus=nf_inq_dimid(ncid,'data_variables',dimidtot)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nvartot',istatus)

          istatus=nf_inq_dimid(ncid,'charsPerLevel',dimchrlvl)
          IF (istatus /= nf_noerr)                                        &
            CALL handle_err('nf_inq_dimid charsPerLevel',istatus)

          istatus=nf_inq_dimid(ncid,'namelen',dimnamlen)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid namelen',istatus)

          istatus=nf_inq_dimid(ncid,'x',dimidnx)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nx',istatus)

          istatus=nf_inq_dimid(ncid,'y',dimidny)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid ny',istatus)

          istatus=nf_inq_dimid(ncid,'levels_1',dimid1)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid 1',istatus)

          istatus=nf_inq_dimid(ncid,nzstr,dimidnz)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid levels',istatus)

          istatus=nf_inq_dimid(ncid,'nav',dimidnav)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_dimid nav',istatus)
!
!  Get variable ids
!
          DO ivar=1,n3dvar
            istatus=nf_inq_varid(ncid,v3dnam(ivar),v3did(ivar))
            IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 3d',istatus)
            ls=LEN(v3dnam(ivar))
            CALL strlnth(v3dnam(ivar),ls)
            WRITE(v3dinv,'(a,a)') v3dnam(ivar)(1:ls),'Inventory'
            istatus=nf_inq_varid(ncid,v3dinv,v3dinvid(ivar))
            IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 3dinv',istatus)
          END DO
          DO ivar=1,n2dvar
            istatus=nf_inq_varid(ncid,v2dnam(ivar),v2did(ivar))
            IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 2d',istatus)
            ls=LEN(v2dnam(ivar))
            CALL strlnth(v2dnam(ivar),ls)
            WRITE(v2dinv,'(a,a)') v2dnam(ivar)(1:ls),'Inventory'
            istatus=nf_inq_varid(ncid,v2dinv,v2dinvid(ivar))
            IF (istatus /= nf_noerr) CALL handle_err('nf_inq_varid 2dinv',istatus)
          END DO

          istatus=nf_inq_varid(ncid,'valtimeMINUSreftime',vtmrefid)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus)

! we're going to have trouble here need to replace with a read of the 
! number of the number of records in the dimensions

          istatus=nf_get_vara_int(ncid,vtmrefid,1,ntwrt,vtmref)
          IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_int',istatus)

          istatus=nf_inq_varid(ncid,'reftime',reftimid)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus)
          istatus=nf_get_vara_double(ncid,reftimid,1,1,reftime)
          IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_double',istatus)

          istatus=nf_inq_varid(ncid,'valtime',valtimid)
          IF (istatus /= nf_noerr) CALL handle_err('nf_inq_var_id',istatus)
          istatus=nf_get_vara_double(ncid,valtimid,1,ntwrt,valtime)
          IF (istatus /= nf_noerr) CALL handle_err('nf_get_vara_double',istatus)
 
        END IF
!
!-----------------------------------------------------------------------
!
!  Set inventory arrays.
!  inventory(k:k) appends the right number of 1's to the string
!
!-----------------------------------------------------------------------
!
        DO k=1,nprlvl
          inventory(k:k)='1'
          IF(iprlvl(k) > 999) THEN
            WRITE(v3dlvl(k),'(a,i4)') 'MB ',iprlvl(k)
          ELSE IF(iprlvl(k) > 99) THEN
            WRITE(v3dlvl(k),'(a,i3)') 'MB ',iprlvl(k)
          ELSE
            WRITE(v3dlvl(k),'(a,i2)') 'MB ',iprlvl(k)
          END IF
        END DO

        ! add surface to 3d level definition (tdo)
        k=ncdflvls
        inventory(k:k)='1'
        WRITE(v3dlvl(k),'(a)') 'SFC'

        !
        !       added a blank inventory (tdo)
        !
        inventory_sfc='1'
        inventory_blank=''

        CALL ctim2abss(year,month,day,hour,minute,second,itimref)
        itimref=itimref-itim1970
!
!-----------------------------------------------------------------------
!
!    End of first-read intializations
!
!-----------------------------------------------------------------------
!
        init=.true.
      END IF
!
!-----------------------------------------------------------------------
!
!    Begin ARPS data conversions
!
!-----------------------------------------------------------------------
!
      itime=itimref+nint(time)
      idxtime=ifile+ntwrt
      valtime(idxtime)=FLOAT(itime)
      vtmref(idxtime)=nint(time)
      reftime(idxtime)=FLOAT(itimref)

      DO i=1,nx
        DO j=1,ny
          rain(i,j)=raing(i,j)+rainc(i,j)
        END DO
      END DO

      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
!
!-----------------------------------------------------------------------
!
      uprt=0.
      vprt=0.
      wprt=0.
      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 and -ln(p) into tem3
!
!-----------------------------------------------------------------------
!
      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

      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
!
!-----------------------------------------------------------------------
!
!  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)
!
!  Store k=2 theta-e and dewpt in tem1,
!  level 1 and 2 respectively.
!
      DO j=1,ny-1
        DO i=1,nx-1
          prmb=0.01*pprt(i,j,2)
          tdew=wmr2td(prmb,(1000.*qvprt(i,j,2)))
          tem1(i,j,1)=oe((tem2(i,j,2)-273.15),tdew,prmb) + 273.15
          tem1(i,j,2)=tdew+273.15
          tem1(i,j,3)=prmb
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!  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.  This is approximately 10m in
!  our configuration
!
!-----------------------------------------------------------------------
!

      ! addition of near surface gph (tdo)
      PRINT *, ' Storing near-sfc gh'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, zps(1,1,2),                      &
                    cdf3d(1,1,ncdflvls,idxgh) )

      ! addition of near surface rh (tdo)
      PRINT *, ' Storing near-sfc rh'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, rh(1,1,2),                       &
                    cdf3d(1,1,ncdflvls,idxrh) )

      PRINT *, ' Storing near-sfc pressure'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,3),                     &
                    cdf2d(1,1,1,idxp) )

      PRINT *,' Storing grid-scale rainfall'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, raing,                           &
                    cdf2d(1,1,1,idxlgsp) )

      PRINT *,' Storing convective rainfall'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, rainc,                           &
                    cdf2d(1,1,1,idxcp) )

      PRINT *,' Storing total accumulated rainfall'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, rain,                            &
                    cdf2d(1,1,1,idxtp) )

     ! delete sfu from netcdf file
     !PRINT *, ' Storing near-sfc u wind'
     !CALL cdf_fill(nx,ny,nxcdf,nycdf, uprt(1,1,2),                     &
     !                    cdf2d(1,1,1,idxsfu) )

      ! add write of u wind into 3d array (tdo)
      PRINT *, ' Storing near-sfc u wind into 3 d array'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, uprt(1,1,2),                     &
                    cdf3d(1,1,ncdflvls,idxuw) )

!      delete sfv from netcdf file
!      PRINT *, ' Storing near-sfc v wind'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, vprt(1,1,2),                     &
!                    cdf2d(1,1,1,idxsfv) )

!     add write of v wind into 3d array (tdo)
      PRINT *, ' Storing near-sfc v wind into 3 d array'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, vprt(1,1,2),                     &
                    cdf3d(1,1,ncdflvls,idxvw) )

!     add write of w wind into 3d array (tdo)
      PRINT *, ' Storing near-sfc w wind into 3 d array'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, wprt(1,1,2),                     &
                    cdf3d(1,1,ncdflvls,idxww) )

!      delete sft from netcdf file
!      PRINT *, ' Storing near-sfc temperature'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2),                     &
!                    cdf2d(1,1,1,idxsft) )

!     add the storage of surface data into the 3d array (tdo)
      PRINT *, ' Storing near-sfc temperature into 3d array'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2),                     &
                      cdf3d(1,1,ncdflvls,idxt) )

      PRINT *, ' Storing near-sfc dew point temperature'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,2),                     &
                    cdf2d(1,1,1,idxdpt) )

      PRINT *,' Storing near-sfc specific humidity'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, qvprt(1,1,2),                    &
                    cdf2d(1,1,1,idxsh) )

      PRINT *, ' Storing near-sfc theta-e'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                     &
                    cdf2d(1,1,1,idxept) )

      PRINT *, ' Storing near-sfc LI'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, li,                              &
                    cdf2d(1,1,1,idxsli) )

      PRINT *, ' Storing near-sfc CAPE'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, pbe,                             &
                    cdf2d(1,1,1,idxcape) )

      PRINT *, ' Storing near-sfc CIN'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, nbe,                             &
                    cdf2d(1,1,1,idxcin) )

!      PRINT *, ' Storing near-sfc LCL'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, lcl,                             &
!                    cdf2d(1,1,1,idxlcl) )

!      PRINT *, ' Storing near-sfc LFC'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, lfc,                             &
!                    cdf2d(1,1,1,idxlfc) )

!      PRINT *, ' Storing near-sfc EL'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, el,                              &
!                    cdf2d(1,1,1,idxel) )

!      PRINT *, ' Storing wet bulb temp diff'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, twdf,                            &
!                    cdf2d(1,1,1,idxtwdf) )

!      PRINT *, ' Storing Cap Strength'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, tcap,                            &
!                    cdf2d(1,1,1,idxtcap) )

!      PRINT *, ' Storing 3-7km Shear'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, shr37,                           &
!                    cdf2d(1,1,1,idxshr37) )

      PRINT *, ' Storing u-storm'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, ustrm,                           &
                    cdf2d(1,1,1,idxustrm) )

      PRINT *, ' Storing v-storm'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, vstrm,                           &
                    cdf2d(1,1,1,idxvstrm) )

!      PRINT *, ' Storing low-level SR flow'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, srlfl,                           &
!                    cdf2d(1,1,1,idxsrlfl) )

!      PRINT *, ' Storing mid-level SR flow'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, srmfl,                           &
!                    cdf2d(1,1,1,idxsrmfl) )

      PRINT *, ' Storing helicity'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, heli,                            &
                    cdf2d(1,1,1,idxheli) )

!      PRINT *, ' Storing Bulk Richardson Number'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, brn,                             &
!                    cdf2d(1,1,1,idxbrn) )

!      PRINT *, ' Storing BRN Shear'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, brnu,                            &
!                    cdf2d(1,1,1,idxbrnu) )

!      PRINT *, ' Storing Boundary Layer Convergence'
!      CALL cdf_fill(nx,ny,nxcdf,nycdf, blcon,                           &
!                    cdf2d(1,1,1,idxblcon) )
!

! J.Case, ENSOC Inc. (9/29/2004) -- Added 2D radar variables to output.

      PRINT *, ' Storing Column Max Reflectivity'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, radar_comp,                      &
                    cdf2d(1,1,1,idxcxr) )

      PRINT *, ' Storing Maximum Radar Echo Tops'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, echo_top,                        &
                    cdf2d(1,1,1,idxmret) )

!-----------------------------------------------------------------------
!
!  Calculate constant pressure level data
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Put temperature into tem2 and -ln(p) into tem3.
!
!-----------------------------------------------------------------------
!
      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
!      ex2=5.2558774
!      DO 355 j=1,ny-1
!      DO 355 i=1,nx-1
!        p00 = 0.01*(pprt(i,j,2))
!        tem1(i,j,1)=p00*
!    :            (1.0+gamma*zps(i,j,2)/tem1(i,j,1))**ex2
! 355     CONTINUE

      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 '
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1,                            &
                    cdf2d(1,1,1,idxmslp) )
!
!-----------------------------------------------------------------------
!
!  Store terrain data.
!
!-----------------------------------------------------------------------
!
      PRINT *, ' Storing terrain data.'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, zp(1,1,2),                       &
                    static2d(1,1,idxtop) )
      PRINT *, ' Storing Coriolis data'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, coriolis,                        &
                    static2d(1,1,idxcor) )
      PRINT *, ' Storing spacing data.'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, spacing,                         &
                    static2d(1,1,idxspa) )
!
!-----------------------------------------------------------------------
!
!    Calculate atmospheric state variables.
!
!-----------------------------------------------------------------------
!
      DO klev=1,nprlvl
        rlnpcdf=-ALOG(100.*FLOAT(iprlvl(klev)))
        PRINT *, ' Storing netCDF data at pr= ',iprlvl(klev)
!
!        Store gh
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     zps,tem3,rlnpcdf,tem1)
        CALL extrph(nx,ny,nz,zps,tem2,pprt,                             &
                    iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1,                          &
                      cdf3d(1,1,klev,idxgh) )
!
!        Store temperature (tem2)
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                      tem2,tem3,rlnpcdf,tem1)
        CALL extrpt(nx,ny,nz,tem2,pprt,zps,                             &
                    iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1,                          &
                      cdf3d(1,1,klev,idxt) )
!
!        Store relative humidity
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     rh,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,rh,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1,                          &
                      cdf3d(1,1,klev,idxrh) )
!
!        Store u and v wind components
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     uprt,tem3,rlnpcdf,tem1(1,1,1))
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     vprt,tem3,rlnpcdf,tem1(1,1,2))
        CALL extrpuv(nx,ny,nz,uprt,vprt,pprt,zps,                       &
                     iprlvl(klev),tem1(1,1,1),tem1(1,1,2))
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxuw) )
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,2),                   &
                      cdf3d(1,1,klev,idxvw) )
!
!        Store vertical velocity
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     wprt,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,wprt,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxww) )

!
!        Store cloud water
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     qc,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,qc,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxcw) )
!
!        Store cloud ice
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     qi,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,qi,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxci) )
!
!        Store turbulent kinetic energy
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     tke,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,tke,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxtk) )
!
! J.Case, ENSCO Inc. (9/29/2004)
!        Store 3D radar reflectivity
!
        CALL v2dinta(nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,                     &
                     radar_3d,tem3,rlnpcdf,tem1)
        CALL extrpq(nx,ny,nz,radar_3d,pprt,iprlvl(klev),tem1)
        CALL cdf_fill(nx,ny,nxcdf,nycdf, tem1(1,1,1),                   &
                      cdf3d(1,1,klev,idxrr) )
      END DO
    END IF   ! Good return from read

    invlen(1)=ncdflvls
    invlen(2)=1
    vecistr(2)=idxtime
    volistart(4)=idxtime
    DO ivar=1,n3dvar
      WRITE(6,'(a,i4,2x,a)') '  Writing 3d ivar:',ivar,                 &
                              v3dlngnam(ivar)

      IF(idxtime == 1) THEN
        istatus=nf_put_vara_text(ncid,v3dlvlid(ivar),                   &
                                 planeistr,lvllen,v3dlvl)
        IF (istatus /= nf_noerr)                                        &
          CALL handle_err('nf_put_vara_text 3d lvl',istatus)
      END IF

      istatus=nf_put_vara_text(ncid,v3dinvid(ivar),                     &
                               vecistr,invlen,inventory)
      IF (istatus /= nf_noerr)                                          &
        CALL handle_err('nf_put_vara_text 3d inv',istatus)

      istatus=nf_put_vara_real(ncid,v3did(ivar),                        &
                   volistart,vollen,cdf3d(1,1,1,ivar))
      IF (istatus /= nf_noerr)                                          &
          CALL handle_err('nf_put_vara_real 3d var',istatus)
    END DO
!
    invlen(1)=1
    invlen(2)=1
    planeistr(4)=idxtime
    DO ivar=1,n2dvar
      WRITE(6,'(a,i4,2x,a)') '  Writing 2d ivar:',ivar,                 &
                                v2dlngnam(ivar)
      IF(adddefs) THEN
        IF(idxtime == 1) THEN
          IF(ivar == idxmslp) THEN
            v2dlvl='MSL'
          ELSE
            v2dlvl='SFC'
          END IF
          istatus=nf_put_vara_text(ncid,v2dlvlid(ivar),                 &
                                 vecistr,sfclen,v2dlvl)
          IF (istatus /= nf_noerr)                                      &
              CALL handle_err('nf_put_vara_text 2d lvl',istatus)
        END IF
      END IF

      istatus=nf_put_vara_text(ncid,v2dinvid(ivar),                     &
                                 vecistr,invlen,inventory_sfc)
      IF (istatus /= nf_noerr)                                          &
            CALL handle_err('nf_put_vara_text 2d inv',istatus)

      istatus=nf_put_vara_real(ncid,v2did(ivar),                        &
                     planeistr,planelen,cdf2d(1,1,1,ivar))
      IF (istatus /= nf_noerr) CALL handle_err('nf_put_vara_real 2d',istatus)
    END DO

  END DO

  IF(adddefs) THEN
!
!  Write static variables
!
    DO ivar=1,nstvar
      WRITE(6,'(a,i4,2x,a)') '  Writing st ivar:',ivar,                 &
                                vstlngnam(ivar)
      istatus=nf_put_vara_real(ncid,vstid(ivar),                        &
                       planeistr,planelen,static2d(1,1,ivar))
      IF (istatus /= nf_noerr)  CALL handle_err('nf_put_vara_real static',istatus)
    END DO

    ls=LEN(origin)
    CALL strlnth(origin,ls)
    istatus=nf_put_vara_text(ncid,origid,1,ls,origin)
    ls=LEN(model)
    CALL strlnth(model,ls)
    istatus=nf_put_vara_text(ncid,modelid,1,ls,model)

!
!   Write 1-d variables if idxtime is 1 (i.e. the first history file)
!
    IF (idxtime == 1) THEN

!
!     set up the first element of the time arrays
!     if more than 1 more valid time will be written to the file
!       foreach time
!         set valtime, vtmref, and reftime using ncdtinterval
!         write out the inventory arrays with blanks
!         AWIPS preprocessor must handle updating inventory arrays if appending
!
      valtime(idxtime)=FLOAT(itime)
      vtmref(idxtime)=nint(time)
      reftime(idxtime)=FLOAT(itimref)

      IF (ncdntime>idxtime) THEN
        DO ntime=idxtime+1,ncdntime
          valtime(ntime)=valtime(idxtime)+FLOAT((ntime-1)*ncdtinterval)
          vtmref(ntime)=vtmref(idxtime)+(ntime-1)*ncdtinterval
          reftime(ntime)=reftime(idxtime)
!
!         write the 3d var inventory with the fill character
!         start array: vecistr(1)=1 vecistr(2)=ntime
!         count array: invlen(1)=1 (1 character) invlen(2)=1
!
          vecistr(1)=1 
          vecistr(2)=ntime
          invlen(1)=1
          invlen(2)=1
          DO ivar=1,n3dvar
            istatus=nf_put_vara_text(ncid,v3dinvid(ivar),                     &
                                  vecistr,invlen,nf_fill_char)
            IF (istatus /= nf_noerr)                                          &
              CALL handle_err('nf_put_vara_text 3d inv',istatus)
          END DO
!
!         write the 2d var inventory with the fill character
!         start array: vecistr(1)=1 vecistr(2)=ntime
!         count array: invlen(1)=1 invlen(2)=1
          DO ivar=1,n2dvar
            istatus=nf_put_vara_text(ncid,v2dinvid(ivar),                     &
                                  vecistr,invlen,nf_fill_char)
            IF (istatus /= nf_noerr)                                          &
              CALL handle_err('nf_put_vara_text 2d inv',istatus)
          END DO
!
        END DO
      ENDIF


      istatus=nf_put_vara_int(ncid,vtmrefid,                                &
                              1,ncdntime,vtmref)
!      istatus=nf_put_vara_double(ncid,reftimid,                             &
!                               1,ncdntime,reftime)
      istatus=nf_put_vara_double(ncid,reftimid,                             &
                               1,1,reftime)
      istatus=nf_put_vara_double(ncid,valtimid,                             &
                               1,ncdntime,valtime)

!    ELSE
!      WRITE(6,'(a,i4)') '  Writing st ivar:',ivar,                     &
    END IF

  END IF
! end idxtime=1

  ntwrt=ntwrt+nhisfile

!
!  Close file
!
  istatus=nf_close(ncid)
  IF (istatus /= nf_noerr) CALL handle_err('nf_close',istatus)
  STOP

  900 CONTINUE
  WRITE(6,'(a)') ' Error reading input data'
  950 CONTINUE
  WRITE(6,'(a)') ' Error setting up netCDF file'
  STOP
END PROGRAM arps2ncdf

!
!##################################################################
!##################################################################
!######                                                      ######
!######                  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, INTENT(IN)  :: nx,ny,nz
  REAL,    INTENT(IN)  :: zps(nx,ny,nz)
  REAL,    INTENT(IN)  :: t(nx,ny,nz)
  REAL,    INTENT(IN)  :: pr(nx,ny,nz)
  INTEGER, INTENT(IN)  :: iprlvl
  REAL,    INTENT(OUT) :: hgtcdf(nx,ny)

  INCLUDE 'phycst.inc'

  REAL, PARAMETER :: gamma = 0.0065     ! degrees/m  lapse rate
  REAL, PARAMETER :: 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                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE extrpt(nx,ny,nz,t,pr,zps,iprlvl,tcdf) 5
!
!-----------------------------------------------------------------------
!
!  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, INTENT(IN)  :: nx,ny,nz
  REAL,    INTENT(IN)  :: t(nx,ny,nz)
  REAL,    INTENT(IN)  :: pr(nx,ny,nz)
  REAL,    INTENT(IN)  :: zps(nx,ny,nz)
  INTEGER, INTENT(IN)  :: iprlvl
  REAL,    INTENT(OUT) :: 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, INTENT(IN)  :: nx,ny,nz
  REAL,    INTENT(IN)  :: q(nx,ny,nz)
  REAL,    INTENT(IN)  :: pr(nx,ny,nz)
  INTEGER, INTENT(IN)  :: iprlvl
  REAL,    INTENT(OUT) :: 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,zps,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, INTENT(IN)  :: nx,ny,nz
  REAL,    INTENT(IN)  :: us(nx,ny,nz)
  REAL,    INTENT(IN)  :: vs(nx,ny,nz)
  REAL,    INTENT(IN)  :: pr(nx,ny,nz)
  REAL,    INTENT(IN)  :: zps(nx,ny,nz)
  INTEGER, INTENT(IN)  :: iprlvl
  REAL,    INTENT(OUT) :: ucdf(nx,ny)
  REAL,    INTENT(OUT) :: 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
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE HANDLE_ERR                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE handle_err(string,istatus) 57
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Report netCDF error in plain English.  Exit.
!
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  CHARACTER(LEN=*), INTENT(IN) :: string
  INTEGER,          INTENT(IN) :: istatus

  CHARACTER (LEN=80) :: nf_strerror
  WRITE(6,'(a,a,/a)') 'netCDF error: ',string, nf_strerror(istatus)

  STOP
END SUBROUTINE handle_err
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE CDF_FILL                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE cdf_fill(nx,ny,nxcdf,nycdf,arpsvar,cdfvar)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Fill netCDF output array which is the physical domain of the
!  ARPS scalar grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: nx,ny
  INTEGER, INTENT(IN)  :: nxcdf,nycdf
  REAL,    INTENT(IN)  :: arpsvar(nx,ny)
  REAL,    INTENT(OUT) :: cdfvar(nxcdf,nycdf)

  INTEGER :: i,j

  DO j=2,ny-2
    DO i=2,nx-2
      cdfvar(i-1,j-1)=arpsvar(i,j)
    END DO
  END DO

  RETURN
END SUBROUTINE cdf_fill