PROGRAM arps2ncdf,70
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  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)
!
!
!-----------------------------------------------------------------------
!
!  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)
!
!    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
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep 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 :: nstyps            ! Maximum number of soil types.

  INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
  PARAMETER (nhisfile_max=200)
  CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'arps2ncdf.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 :: 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 :: tsfc   (:,:,:) ! Temperature at surface (K)
  REAL, ALLOCATABLE :: tsoil  (:,:,:) ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: wetsfc (:,:,:) ! Surface soil moisture
  REAL, ALLOCATABLE :: wetdp  (:,:,:) ! Deep 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 :: 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(:,:,:)

  REAL, ALLOCATABLE :: mf2d(:,:)
  REAL, ALLOCATABLE :: coriolis(:,:)
  REAL, ALLOCATABLE :: spacing(:,:)
!
!-----------------------------------------------------------------------
!
!  Misc ARPS variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nch
  REAL :: time
  REAL :: latnot(2)
  LOGICAL :: init
!
!-----------------------------------------------------------------------
!
!  netCDF dimensions
!
!  nprlvl is difined in arps2cdf.inc
!
!-----------------------------------------------------------------------
!
  INTEGER :: nxcdf,nycdf,mxtime
  PARAMETER(mxtime=40)
  INTEGER :: n3dvar,n2dvar,nvartot,nstvar
  PARAMETER (n3dvar=6,                                                  &
             n2dvar=28,                                                 &
             nvartot=(n3dvar+n2dvar),                                   &
             nstvar=3)
!
!-----------------------------------------------------------------------
!
!  netCDF variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: istatus,ncid
  INTEGER :: dimidrec,dimidtot,dimidnt,dimidnav
  INTEGER :: dimidnx,dimidny,dimid1,dimidnz,dimchrlvl,dimnamlen
  INTEGER :: vtmrefid,valtimid,reftimid,origid,modelid
!
  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)
!
  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)
  INTEGER :: ntime
!
  CHARACTER (LEN=4) :: v3dnam(n3dvar)
  CHARACTER (LEN=4) :: v2dnam(n2dvar)
  CHARACTER (LEN=14) :: vstnam(nstvar)
  CHARACTER (LEN=10) :: v3dlvl(nprlvl)
  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) :: inventory
  CHARACTER (LEN=1) :: inventory_sfc
  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
  PARAMETER (idxgh=1,idxrh=2,idxt=3,                                    &
             idxuw=4,idxvw=5,idxww=6)
!
  DATA v3dnam /'gh','rh','t','uw','vw','ww'/
  DATA v3dlngnam /'Geopotential Height',                                &
                  'Relative Humidity',                                  &
                  'Temperature',                                        &
                  'u Wind Component',                                   &
                  'v Wind Component',                                   &
                  'w Wind Component'/
  DATA v3dunits /'geopotential meters',                                 &
                 'percent',                                             &
                 'degrees K',                                           &
                 'meters/second',                                       &
                 'meters/second',                                       &
                 'meters/second'/
  DATA v3dmin /    0.,  0.,150.,-300.,-300.,-100./
  DATA v3dmax /40000.,110.,400., 300., 300., 100./
!
  INTEGER :: idxmslp,idxp,idxlgsp,idxcp,idxtp
  INTEGER :: idxsfu,idxsfv,idxsft,idxdpt
  INTEGER :: idxsh,idxept,idxsli,idxcape,idxcin
  INTEGER :: idxlcl,idxlfc,idxel,idxtwdf,idxtcap
  INTEGER :: idxshr37,idxustrm,idxvstrm
  INTEGER :: idxsrlfl,idxsrmfl,idxheli,idxbrn,idxbrnu,idxblcon
  PARAMETER (idxmslp=1,                                                 &
             idxp   =2,                                                 &
             idxlgsp=3,                                                 &
             idxcp  =4,                                                 &
             idxtp  =5,                                                 &
             idxsfu =6,                                                 &
             idxsfv =7,                                                 &
             idxsft =8,                                                 &
             idxdpt =9,                                                 &
             idxsh  =10,                                                &
             idxept =11,                                                &
             idxsli =12,                                                &
             idxcape=13,                                                &
             idxcin =14,                                                &
             idxlcl =15,                                                &
             idxlfc =16,                                                &
             idxel  =17,                                                &
             idxtwdf =18,                                               &
             idxtcap =19,                                               &
             idxshr37=20,                                               &
             idxustrm=21,                                               &
             idxvstrm=22,                                               &
             idxsrlfl=23,                                               &
             idxsrmfl=24,                                               &
             idxheli =25,                                               &
             idxbrn  =26,                                               &
             idxbrnu =27,                                               &
             idxblcon=28)

  DATA v2dnam /'mslp','p','lgsp','cp','tp','sfu','sfv',                 &
               'sft','dpt','sh','ept','sli','cape','cin','lcl',         &
               'lfc','el','twdf','tcap','shr37','ustrm','vstrm',        &
               'srlfl','srmfl','heli','brn','brnu','blcon'/

  DATA v2dlngnam /'Mean Sea Level Pressure',                            &
                  'Surface Pressure',                                   &
                  'Grid Scale Precipitation',                           &
                  'Convective Precipitation',                           &
                  'Total Precipitation',                                &
                  'Surface u Wind Component',                           &
                  'Surface v Wind Component',                           &
                  'Surface Temperature',                                &
                  'Surface Dew Point',                                  &
                  'Specific Humidity',                                  &
                  'Theta-e',                                            &
                  'Surface Lifted Index',                               &
                  'CAPE',                                               &
                  'CIN',                                                &
                  'Lifting Condensation Lvl',                           &
                  'Level of Free Convection',                           &
                  'Equilibrium Level',                                  &
                  'Max Wet-Bulb Temp Diff',                             &
                  'Cap Strength',                                       &
                  '3-7km Shear',                                        &
                  'Est Storm Motion u comp',                            &
                  'Est Storm Motion v comp',                            &
                  'Storm Rel Low-Level Flow',                           &
                  'Storm Rel Mid-Level Flow',                           &
                  'Storm Rel Helicity',                                 &
                  'Bulk Richardson Number',                             &
                  'BRN Shear u',                                        &
                  'Bound-Layer Conv x 1000'/

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

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

  DATA v2dmax / 110000.,120000.,1000.,1000.,1000.,                      &
                100.,100.,                                              &
                100.,100.,1.,                                           &
                500.,100.,10000.,10000.,                                &
                40000.,40000.,40000.,100.,100.,1.,100.,100.,            &
                100.,100.,900.,900.,100.,10./
!
  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=20) :: 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 :: itim1970
  PARAMETER (itim1970=315619200)
  CHARACTER (LEN=8) :: cdldate
  CHARACTER (LEN=40) :: depictor
  CHARACTER (LEN=132) :: origin,model
  PARAMETER (cdldate='19990301',                                        &
             depictor='TEST_ARPS',                                      &
             origin='HubCAPS',                                          &
             model='ARPS 4.4.1')
!
!-----------------------------------------------------------------------
!
!  Namelists
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=132) :: cdfname
  INTEGER :: ipktyp,nbits
  NAMELIST /ncdf_options/ cdfname,ipktyp,nbits
!
!-----------------------------------------------------------------------
!
!  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
!
  INTEGER :: i,j,k,klev,ifile,ireturn,imid,jmid
  REAL :: rlnpcdf,r2d
  REAL :: ctrx,ctry,swx,swy
  REAL :: gamma,ex1,ex2,rln700,p00,t00
  REAL :: prmb,tdew

  REAL, ALLOCATABLE :: p_mb(:,:,:)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Initializations
!
!-----------------------------------------------------------------------
!
!
  WRITE(6,'(6(/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.
!
!-----------------------------------------------------------------------
!
  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,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

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

  planelen(1)=nxcdf
  planelen(2)=nycdf
  planelen(3)=1
  planelen(4)=1
  lvllen(1)=10
  lvllen(2)=nprlvl
  sfclen(1)=10
  sfclen(2)=1
  vollen(1)=nxcdf
  vollen(2)=nycdf
  vollen(3)=nprlvl
  vollen(4)=1
!
!-----------------------------------------------------------------------
!
!  Get output options and the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!
  ipktyp=1
  nbits=16
  READ(5,ncdf_options,END=900)

  ntime=nhisfile

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

  ALLOCATE(x      (nx))
  ALLOCATE(y      (ny))
  ALLOCATE(z      (nz))
  ALLOCATE(zp     (nx,ny,nz))

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

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

  ALLOCATE(tsfc   (nx,ny,0:nstyps))
  ALLOCATE(tsoil  (nx,ny,0:nstyps))
  ALLOCATE(wetsfc (nx,ny,0:nstyps))
  ALLOCATE(wetdp  (nx,ny,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(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))

  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,nprlvl,n3dvar))

  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

  tsfc   =0.0
  tsoil  =0.0
  wetsfc =0.0
  wetdp  =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

  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

  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.

  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,nstyps,hinfmt,nch,grdbasfn(1:lengbf),         &
               lengbf,                                                  &
               hisfile(ifile)(1:lenfil),lenfil,time,                    &
               x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,                &
               qvprt, qc, qr, qi, qs, qh, tke, kmh, kmv,                &
               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,            &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,                                      &
               usflx,vsflx,ptsflx,qvsflx,                               &
               ireturn, tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
    IF( ireturn == 0 ) THEN   ! successful read
      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


!       nf_noerr = 0 ! correct??
!
!-----------------------------------------------------------------------
!
!  Open netCDF file
!
!-----------------------------------------------------------------------
!
        istatus=nf_create(cdfname,nf_clobber,ncid)
        IF (istatus /= nf_noerr) CALL handle_err('nf_create',istatus)
!
!-----------------------------------------------------------------------
!
!  Define dimensions
!
!-----------------------------------------------------------------------
!
        IF(nprlvl < 10) THEN
          WRITE(nzstr,'(a,i1)') 'levels_',nprlvl
        ELSE IF(nprlvl < 100) THEN
          WRITE(nzstr,'(a,i2)') 'levels_',nprlvl
        ELSE
          WRITE(nzstr,'(a,i3)') 'levels_',nprlvl
        END IF
!
        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',ntime,dimidnt)
        IF (istatus /= nf_noerr) CALL handle_err('nf_def_dim ntime',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,nprlvl,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)
!
        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)=dimidrec
!
!-----------------------------------------------------------------------
!
!  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,nprlvl)
          istatus=nf_put_att_text(ncid,v3did(ivar),'levels',            &
                          2,'MB')
!
          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)
!
          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)=dimidrec
        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)
!
          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,1,                  &
                           dimidnt,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)
!
!  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

!
!-----------------------------------------------------------------------
!
!  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 = ATAN(1.0)*4.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
!
        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)
!
!-----------------------------------------------------------------------
!
!  Set inventory arrays.
!
!-----------------------------------------------------------------------
!
        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

!
        inventory_sfc='1'
!
        CALL ctim2abss(year,month,day,hour,minute,second,itimref)
        itimref=itimref-itim1970
        PRINT *, ' reference time: ',itimref
        DO i=1,ntime
          reftime(i)=FLOAT(itimref)
        END DO
!
!-----------------------------------------------------------------------
!
!    End of first-read intializations
!
!-----------------------------------------------------------------------
!
        init=.true.
      END IF
!
!-----------------------------------------------------------------------
!
!    Begin ARPS data conversions
!
!-----------------------------------------------------------------------
!
      itime=itimref+nint(time)
      valtime(ifile)=FLOAT(itime)
      vtmref(ifile)=nint(time)

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

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

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

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

      PRINT *, ' Storing near-sfc temperature'
      CALL cdf_fill(nx,ny,nxcdf,nycdf, tem2(1,1,2),                     &
                    cdf2d(1,1,1,idxsft) )

      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) )
!
!-----------------------------------------------------------------------
!
!  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) )
      END DO
    END IF   ! Good return from read
!
    invlen(1)=nprlvl
    invlen(2)=1
    vecistr(2)=ifile
    volistart(4)=ifile
    DO ivar=1,n3dvar
      WRITE(6,'(a,i4,2x,a)') '  Writing 3d ivar:',ivar,                 &
                              v3dlngnam(ivar)
!
      IF(ifile == 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)=ifile
    DO ivar=1,n2dvar
      WRITE(6,'(a,i4,2x,a)') '  Writing 2d ivar:',ivar,                 &
                              v2dlngnam(ivar)
      IF(ifile == 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

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

!
!  Write 1-d variables
!
  istatus=nf_put_vara_int(ncid,vtmrefid,                                &
                             1,ntime,vtmref)
  istatus=nf_put_vara_double(ncid,reftimid,                             &
                             1,ntime,reftime)
  istatus=nf_put_vara_double(ncid,valtimid,                             &
                             1,ntime,valtime)
  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)
!
!  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(nx,ny,nz,zps,t,pr,iprlvl,hgtcdf) 2
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE EXTRPH                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Extrapolate height by using a std atmos lapse rate
!  below the last physical level.  Above the domain,
!  assume a constant temperature above 300 mb, otherwise
!  use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: zps(nx,ny,nz)
  REAL :: t(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: hgtcdf(nx,ny)
!
  INCLUDE 'phycst.inc'
!
  REAL :: gamma,rddg,const
  PARAMETER ( gamma = 0.0065,    & ! degrees/m  lapse rate
          rddg  = (rd/g),                                               &
              const = (rd*gamma/g) )
!
  INTEGER :: i,j
  REAL :: prcdf
!
  prcdf=100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf < pr(i,j,nz-1)) THEN
        IF(pr(i,j,nz-1) <= 30000.) THEN
          hgtcdf(i,j)=zps(i,j,nz-1) +                                   &
              rddg*t(i,j,nz-1)*ALOG(pr(i,j,nz-1)/prcdf)
        ELSE
          hgtcdf(i,j)=zps(i,j,nz-1) + (t(i,j,nz-1)/gamma)*              &
                      (1.-(prcdf/pr(i,j,nz-1))**const)
        END IF
      ELSE IF(prcdf >= pr(i,j,2)) THEN
        hgtcdf(i,j)=zps(i,j,2) + (t(i,j,2)/gamma)*                      &
               (1.-(prcdf/pr(i,j,2))**const)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrph
!


SUBROUTINE extrpt(nx,ny,nz,t,pr,zps,iprlvl,tcdf) 3
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE EXTRPT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Extrapolate temperature by using a std atmos lapse rate
!  below the last physical level.  Above the domain,
!  assume a constant temperature above 300 mb, otherwise
!  use the std atmos lapse rate.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  REAL :: t(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  REAL :: zps(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: tcdf(nx,ny)
!
  INCLUDE 'phycst.inc'
!
  REAL :: gamma,const
  PARAMETER ( gamma = 0.0065,    & ! degrees/m  lapse rate
          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(nx,ny,nz,q,pr,iprlvl,qcdf) 8
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE EXTRPQ                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign a zero-gradient value to scalars below ground.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster, from arps2gem
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: q(nx,ny,nz)
  REAL :: pr(nx,ny,nz)
  INTEGER :: iprlvl
  REAL :: qcdf(nx,ny)
!
  INTEGER :: i,j
  REAL :: prcdf
!
  prcdf=100.*FLOAT(iprlvl)
  DO j=1,ny-1
    DO i=1,nx-1
      IF(prcdf < pr(i,j,nz-1)) THEN
        qcdf(i,j)=q(i,j,nz-1)
      ELSE IF(prcdf > pr(i,j,2)) THEN
        qcdf(i,j)=q(i,j,2)
      END IF
    END DO
  END DO
  RETURN
END SUBROUTINE extrpq


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


SUBROUTINE handle_err(string,istatus) 32
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE HANDLE_ERR                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Report netCDF error in plain English.  Exit.
!
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster
!  2/19/99
!
!  MODIFICATION HISTORY:
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  CHARACTER (LEN=*) :: string
  CHARACTER (LEN=80) :: nf_strerror
  INTEGER :: istatus
  WRITE(6,'(a,a,/a)') 'netCDF error: ',string,                          &
                      nf_strerror(istatus)
  STOP
END SUBROUTINE handle_err
!


SUBROUTINE cdf_fill(nx,ny,nxcdf,nycdf,                                  &
           arpsvar,cdfvar)
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE CDF_FILL                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  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 :: nx,ny
  INTEGER :: nxcdf,nycdf
  REAL :: arpsvar(nx,ny)
  REAL :: 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