PROGRAM arpsextsnd,86
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM EXTSND                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Extracts soundings from ARPS model history data.
!  User specifies an ARPS history file and either a lat,lon or
!  an x-y location.
!
!  Creates a file with the sounding data in an ASCII format.
!  The format matches one used by UNIDATA's GEMPAK software for ASCII
!  output of sounding data (GEMPAK program SNLIST).
!
!  Because of that compatibility, the output file can be
!  read by certain sounding plotting and analysis programs on the
!  computers at the University of Oklahoma (OU).
!
!  The following instructions apply to ARPS users at OU:
!
!  After running this program use
!  ~rcarpent/bin/skewt -sfc -hodo data_file
!
!  where data_file is the output of this program (a name specified
!  in the input of this program).  An NCARgraphics gmeta file
!  is created by skewt.
!
!  It is also possible to use this file as input to the standard
!  skewt plotting programs on the metgem and Rossby computers.
!  See the author for details.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  April 1992.
!
!  MODIFICATION HISTORY:
!
!  October, 1992 (K. Brewster)
!  Conversion to ARPS 3.0.
!
!  February, 1993 (K. Brewster)
!  Additional documentation for ARPS 3.1 release
!
!  4/13/93 (M. Xue)
!  Modified to conform to the new data dump format.
!
!  10/19/94 (KB)
!  Modified to conform to yet another data dump format.
!  Corrections made for change in x,y definition in dump file.
!
!  06/19/95 (KB)
!  Modified documentation and output formats to update info.
!
!  03/07/96 (KB)
!  Changed calling arguments in FINDLC, moved subroutines to a library.
!  Corrected a bug in ymap(ny) calculation.
!
!  11/06/96 (KB)
!  Unified interpolation calls with other programs using interpolation.
!  Added capability to process multiple soundings in one run.
!
!  2000/05/19 (Gene Bassett)
!    Converted to F90, creating allocation and arpsextsnd main
!    subroutines.
!
!  2000/07/28 (Ming Xue)
!  Converted to F90 free format. Use ALLOCATABLE instead of
!  POINTER allocation to avoid double memory usage.
!
!  Changed to namelist format input.
!
!  2001/12/05 (Ming Hu)
!    Add an output file which is sounding used in
!         ARPS as base field like May20.snd
!
!-----------------------------------------------------------------------
!
!  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)
!
!    uprt     Perturbation x component of velocity (m/s)
!    vprt     Perturbation y component of velocity (m/s)
!    wprt     Perturbation z component of velocity (m/s)
!
!    ptprt    Perturbation potential temperature (K)
!    pprt     Perturbation pressure (Pascal)
!
!    qvprt    Perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state 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
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  CALCULATED DATA ARRAYS:
!
!    su       Sounding x component of velocity (m/s)
!    sv       Sounding y component of velocity (m/s)
!    stheta   Sounding potential temperature (K)
!    spres    Sounding pressure (mb)
!    stemp    Sounding temperature (C)
!    sdewp    Sounding dew-point (C)
!    sdrct    Sounding wind direction (degrees)
!    ssped    Sounding wind speed (m/s)
!    shght    Sounding height (m)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!
!   Temporary arrays are defined and used differently by each
!   subroutine.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz

  INTEGER :: nstyps
  INTEGER, PARAMETER :: maxsnd = 200
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.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 on the staggered grid.
  REAL, ALLOCATABLE :: zp    (:,:,:)     ! The physical height coordinate defined at
                                         ! w-point of the staggered grid.

  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 :: 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
                                         ! humidity (kg/kg)

  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 :: 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)
!
!-----------------------------------------------------------------------
!
!  Computed variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: xs(:)      ! x location of scalar points
  REAL, ALLOCATABLE :: ys(:)      ! y location of scalar points
!
!-----------------------------------------------------------------------
!
!  Extracted sounding variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: su(:,:)
  REAL, ALLOCATABLE :: sv(:,:)
  REAL, ALLOCATABLE :: stheta(:,:)
  REAL, ALLOCATABLE :: sqv(:,:)
  REAL, ALLOCATABLE :: spres(:,:)
  REAL, ALLOCATABLE :: stemp(:,:)
  REAL, ALLOCATABLE :: sdewp(:,:)
  REAL, ALLOCATABLE :: sdrct(:,:)
  REAL, ALLOCATABLE :: ssped(:,:)
  REAL, ALLOCATABLE :: shght(:,:)
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: dxfld(:)
  REAL, ALLOCATABLE :: dyfld(:)
  REAL, ALLOCATABLE :: rdxfld(:)
  REAL, ALLOCATABLE :: rdyfld(:)
  REAL, ALLOCATABLE :: slopey(:,:,:)
  REAL, ALLOCATABLE :: alphay(:,:,:)
  REAL, ALLOCATABLE :: betay(:,:,:)
  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)

  REAL :: xpt(maxsnd),ypt(maxsnd)
  CHARACTER (LEN=256) :: ofile(maxsnd)
  INTEGER :: ipt(maxsnd),jpt(maxsnd)
!
!-----------------------------------------------------------------------
!
!  Sounding indentification variables
!  These are required for the proper operation of the
!  plotting programs on the metgem computer at the
!  University of Oklahoma.
!
!-----------------------------------------------------------------------
!
!wdt update
  CHARACTER (LEN=8) :: stid(maxsnd)
  REAL :: slat(maxsnd),slon(maxsnd),selev(maxsnd)
  INTEGER :: istnm(maxsnd)
!
!-----------------------------------------------------------------------
!
!  Map projection variables
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: xmap(:)
  REAL, ALLOCATABLE :: ymap(:)
  REAL, ALLOCATABLE :: latgr(:,:)
  REAL, ALLOCATABLE :: longr(:,:)

  INTEGER :: i4time,iyr,imo,iday,ihr,imin,isec
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc. internal variables
!
!-----------------------------------------------------------------------
!
  REAL :: time,xctr,yctr,xll,yll,xsndmap,ysndmap
  REAL :: ustorm,vstorm
  INTEGER :: i,j,k,locopt,isnd,nsnd
  INTEGER :: ireturn,hinfmt,lengbf,lenfil,nchin
  CHARACTER (LEN=256) :: filename,grdbasfn
  CHARACTER (LEN=80) :: sndtime,snddate,sndlocation

  NAMELIST /input/ hinfmt,grdbasfn,filename,ustorm,vstorm,              &
                  locopt,nsnd,slat,slon,xpt,ypt,ofile,stid,istnm

  INTEGER :: nlevs, istatus
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  WRITE(6,'(/6(/1x,a)/)')                                               &
      '###############################################################',&
      '#                                                             #',&
      '# Welcome to EXTSND, a program that reads in a history file   #',&
      '# generated by ARPS and extracts a sounding at an x-y point.  #',&
      '#                                                             #',&
      '###############################################################'

  xpt = -999.0
  ypt = -999.0

  READ(5,input)

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

  lenfil = LEN_trim(filename)
  WRITE(6,'(a,a)')' The data set name is ', filename(1:lenfil)

  IF( nsnd > maxsnd )  then
    WRITE(6,'(a,/a,i5)')                                                &
    'The number of sounding locations to be extracted exceeded maximum ',&
    'allowed. nsnd is reset to ', maxsnd
    nsnd = maxsnd
  ENDIF

!
!-----------------------------------------------------------------------
!
!  Obtain the grid dimensions from input data.
!
!-----------------------------------------------------------------------
!

  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
!
!-----------------------------------------------------------------------
!
! Allocate arrays  
!
!-----------------------------------------------------------------------
!
  ALLOCATE(x      (nx),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:x')
  ALLOCATE(y      (ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:y')
  ALLOCATE(z      (nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:z')
  ALLOCATE(zp     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:zp')
  ALLOCATE(uprt   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:uprt')
  ALLOCATE(vprt   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:vprt')
  ALLOCATE(wprt   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:wprt')
  ALLOCATE(ptprt  (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ptprt')
  ALLOCATE(pprt   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:pprt')
  ALLOCATE(qvprt  (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qvprt')
  ALLOCATE(qc     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qc')
  ALLOCATE(qr     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qr')
  ALLOCATE(qi     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qi')
  ALLOCATE(qs     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qs')
  ALLOCATE(qh     (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qh')
  ALLOCATE(tke    (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tke')
  ALLOCATE(kmh    (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:kmh')
  ALLOCATE(kmv    (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:kmv')
  ALLOCATE(ubar   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ubar')
  ALLOCATE(vbar   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:vbar')
  ALLOCATE(wbar   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:wbar')
  ALLOCATE(ptbar  (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ptbar')
  ALLOCATE(pbar   (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:pbar')
  ALLOCATE(rhobar (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:rhobar')
  ALLOCATE(qvbar  (nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qvbar')

  ALLOCATE(soiltyp (nx,ny,nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:soiltyp')
  ALLOCATE(stypfrct(nx,ny,nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:stypfrct')
  ALLOCATE(vegtyp(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:vegtyp')
  ALLOCATE(lai    (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:lai')
  ALLOCATE(roufns (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:roufns')
  ALLOCATE(veg    (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:veg')
  ALLOCATE(tsfc   (nx,ny,0:nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tsfc')
  ALLOCATE(tsoil  (nx,ny,0:nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tsoil')
  ALLOCATE(wetsfc (nx,ny,0:nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:wetsfc')
  ALLOCATE(wetdp  (nx,ny,0:nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:wetdp')
  ALLOCATE(wetcanp(nx,ny,0:nstyps),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:wetcanp')
  ALLOCATE(snowdpth(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:snowdpth')
  ALLOCATE(raing(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:raing')
  ALLOCATE(rainc(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:rainc')
  ALLOCATE(prcrate(nx,ny,4),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:prcrate')
  ALLOCATE(radfrc(nx,ny,nz),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:radfrc')
  ALLOCATE(radsw (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:radsw')
  ALLOCATE(rnflx (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:rnflx')
  ALLOCATE(usflx (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:usflx')
  ALLOCATE(vsflx (nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:vsflx')
  ALLOCATE(ptsflx(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ptsflx')
  ALLOCATE(qvsflx(nx,ny),stat=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:qvsflx')

  ALLOCATE(xs(nx),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:xs')
  ALLOCATE(ys(ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ys')
  ALLOCATE(su(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:su')
  ALLOCATE(sv(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:sv')
  ALLOCATE(stheta(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:stheta')
  ALLOCATE(sqv(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:sqv')
  ALLOCATE(spres(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:spres')
  ALLOCATE(stemp(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:stemp')
  ALLOCATE(sdewp(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:sdewp')
  ALLOCATE(sdrct(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:sdrct')
  ALLOCATE(ssped(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ssped')
  ALLOCATE(shght(nz,maxsnd),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:shght')
  ALLOCATE(dxfld(nx),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:dxfld')
  ALLOCATE(dyfld(ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:dyfld')
  ALLOCATE(rdxfld(nx),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:rdxfld')
  ALLOCATE(rdyfld(ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:rdyfld')
  ALLOCATE(slopey(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:slopey')
  ALLOCATE(alphay(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:alphay')
  ALLOCATE(betay(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:betay')
  ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tem1')
  ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tem2')
  ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:tem3')

  ALLOCATE(xmap(nx),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:xmap')
  ALLOCATE(ymap(ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:ymap')
  ALLOCATE(latgr(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:latgr')
  ALLOCATE(longr(nx,ny),STAT=istatus)
  CALL check_alloc_status(istatus, 'arpsextsnd:longr')

  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
  radsw =0.0
  qvbar  =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
  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

  xs=0.0
  ys=0.0
  su=0.0
  sv=0.0
  stheta=0.0
  sqv=0.0
  spres=0.0
  stemp=0.0
  sdewp=0.0
  sdrct=0.0
  ssped=0.0
  shght=0.0
  dxfld=0.0
  dyfld=0.0
  rdxfld=0.0
  rdyfld=0.0
  slopey=0.0
  alphay=0.0
  betay=0.0
  tem1=0.0
  tem2=0.0
  tem3=0.0

!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
  CALL dtaread(nx,ny,nz,nstyps,                                         &
               hinfmt,nchin,grdbasfn(1:lengbf),lengbf,                  &
               filename(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

!
!-----------------------------------------------------------------------
!
!  Calculate scalar locations
!
!-----------------------------------------------------------------------
!
    DO i=1,nx-1
      xs(i)=0.5*(x(i)+x(i+1))
    END DO
    DO j=1,ny-1
      ys(j)=0.5*(y(j)+y(j+1))
    END DO
    dx=x(2)-x(1)
    dy=y(2)-y(1)
!
!-----------------------------------------------------------------------
!
!  Compute sounding time and write it
!
!-----------------------------------------------------------------------
!
    CALL ctim2abss( year,month,day,hour,minute,second,i4time)
    i4time=i4time + nint(time)
    CALL abss2ctim( i4time, iyr, imo, iday, ihr, imin, isec )
    WRITE(6,'(a,i4.4,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2,a1,i2.2)')         &
        '  Time of history data: ',                                     &
        iyr,'/',imo,'/',iday,':',ihr,':',imin,':',isec
!
!-----------------------------------------------------------------------
!
!  Find and write the lot, lon of the extraction point in
!  SNLIST file.
!
!-----------------------------------------------------------------------
!
    latnot(1)=trulat1
    latnot(2)=trulat2
    CALL setmapr(mapproj,sclfct,latnot,trulon)
    CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
    PRINT *, ' dx= ',dx,' dy= ',dy
    xll=xctr-(0.5*(nx-3)*dx)
    yll=yctr-(0.5*(ny-3)*dy)

    DO i=1,nx-1
      xmap(i)=xll+xs(i)
    END DO
    xmap(nx)=2.*xmap(nx-1)-xmap(nx-2)
    DO j=1,ny-1
      ymap(j)=yll+ys(j)
    END DO
    ymap(ny)=2.*ymap(ny-1)-ymap(ny-2)
    CALL xytoll(nx,ny,xmap,ymap,latgr,longr)


    DO isnd=1,nsnd
      IF(locopt == 1) THEN

        IF(slat(isnd) < -90.) EXIT

        CALL lltoxy(1,1,slat(isnd),slon(isnd),xpt(isnd),ypt(isnd))
        xpt(isnd)=xpt(isnd)-xll
        ypt(isnd)=ypt(isnd)-yll

      ELSE
        xpt(isnd)=xpt(isnd)*1000.
        ypt(isnd)=ypt(isnd)*1000.
        xsndmap=xpt(isnd)+xll
        ysndmap=ypt(isnd)+yll
        CALL xytoll(1,1,xsndmap,ysndmap,slat(isnd),slon(isnd))
      END IF

      print*,' isnd, xpt, ypt =', isnd, xpt(isnd), ypt(isnd)
!
!-----------------------------------------------------------------------
!
!  Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
      CALL setijloc(1,1,nx,ny,xpt(isnd),ypt(isnd),                      &
                    xs,ys,ipt(isnd),jpt(isnd))
!
      WRITE (6,'(2x,a,f12.2,a,f12.2,/2X,a,f12.2,a,f12.2,/2X,a,i6,a,i6)')&
          ' location x= ',(0.001*xpt(isnd)),'   y= ',(0.001*ypt(isnd)), &
          '        lat= ',slat(isnd),' lon= ',slon(isnd),               &
          ' found near pt i= ',ipt(isnd),'   j= ',jpt(isnd)
      WRITE (6,'(12x,a,f12.2,a,f12.2,/10X,a,f12.2,a,f12.2)')            &
          ' x= ',(0.001*xs(ipt(isnd))),                                 &
          ' y= ',(0.001*ys(jpt(isnd))),                                 &
          ' lat= ',latgr(ipt(isnd),jpt(isnd)),                          &
          ' lon= ',longr(ipt(isnd),jpt(isnd))

      IF(ofile(isnd)(1:4) == '    ') ofile(isnd)='SNLIST.FIL'

    END DO
!
!-----------------------------------------------------------------------
!
!  Interpolate (in the horizontal) for the whole vertical column.
!
!-----------------------------------------------------------------------
!
    CALL setdxdy(nx,ny,                                                 &
                 1,nx-1,1,ny-1,                                         &
                 xs,ys,dxfld,dyfld,rdxfld,rdyfld)
    CALL colint(nx,ny,nz,maxsnd,nsnd,                                   &
                xs,ys,zp,xpt,ypt,ipt,jpt,                               &
                uprt, vprt, ptprt, pprt, qvprt,                         &
                ubar, vbar, ptbar, pbar, qvbar,                         &
                su,sv,stheta,spres,shght,sqv,selev,                     &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                tem1,nlevs)
!
!-----------------------------------------------------------------------
!
!  Add back storm motion that ARPS subtracts from the sounding winds.
!
!-----------------------------------------------------------------------
!
    DO isnd=1,nsnd
      DO k=1,nlevs
        su(k,isnd)=su(k,isnd)+ustorm
        sv(k,isnd)=sv(k,isnd)+vstorm
      END DO
!
!-----------------------------------------------------------------------
!
!  Convert sounding to desired units/quantities
!
!-----------------------------------------------------------------------
!
      CALL cnvsnd(su(1,isnd),sv(1,isnd),stheta(1,isnd),                 &
                spres(1,isnd),sqv(1,isnd),slon(isnd),                   &
                sdrct(1,isnd),ssped(1,isnd),                            &
                stemp(1,isnd),sdewp(1,isnd),nlevs)
!
!-----------------------------------------------------------------------
!
!  Output sounding to look like GEMPAK SNLIST.FIL
!  to use the Skew-T and Hodograph programs
!  attributed to Bill McCaul.
!
!  Example of what the header looks like:
!
!23456789012345678901234567890123456789012345678901234567890
!SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED
!
!STID = SEP        STNM =    72260   TIME = 920308/1500
!SLAT =    32.21   SLON =   -98.18   SELEV =   399.0
!
!  PRES     HGHT     TMPC     DWPC     DRCT     SPED
!
!-----------------------------------------------------------------------
!
      OPEN(3,FILE=trim(ofile(isnd)),STATUS='unknown')
      WRITE(3,'(/a/)') ' SNPARM = PRES;HGHT;TMPC;DWPC;DRCT;SPED'
      iyr=MOD(iyr,100)
!wdt update
      WRITE(3,'(a,a8,3x,a,i8,3x,a,i2.2,i2.2,i2.2,a1,i2.2,i2.2)')        &
          ' STID = ',stid(isnd),                                        &
          'STNM = ',istnm(isnd),                                        &
          'TIME = ',iyr,imo,iday,'/',ihr,imin
      WRITE(3,'(a,f8.2,3x,a,f8.2,3x,a,f7.1/)')                          &
          ' SLAT = ',slat(isnd),'SLON = ',slon(isnd),                   &
          'SELV = ',selev(isnd)
      WRITE(3,'(6x,a)')                                                 &
          'PRES     HGHT     TMPC     DWPC     DRCT     SPED'
!
!-----------------------------------------------------------------------
!
!    Print out of sounding
!
!-----------------------------------------------------------------------
!
      DO k=1,nlevs
        WRITE(3,'(1x,6(F9.2))')                                         &
              spres(k,isnd),shght(k,isnd),                              &
              stemp(k,isnd),sdewp(k,isnd),                              &
              sdrct(k,isnd),ssped(k,isnd)
      END DO
      CLOSE(3)
      write(6,'(a,a,a)')                                                &
      'Sounding file ',trim(ofile(isnd)),' was produced.'
            
!-----------------------------------------------------------------------
!
!  OUTPUT sounding according to follow format which can be used in 
!         ARPS as base field like May20.snd
!
!              Sounding File Format for ARPS 3.0
!
!  Record 1: a header line, such as "1-D Sounding Input for ARPS"
!            (skipped)
!  Record 2: miscellaneous description of sounding (skipped)
!  Record 3: time of sounding (character*72)
!  Record 4: date of sounding (character*72)
!  Record 5: location of sounding (character*72)
!  Record 6: three character strings designate the sounding
!            data type, e.g.
!            'pressure' 'potential_temperature' 'relative_humidity'.
!            Only the first character of the strings
!            is decoded, and thus the first character should not be
!            left blank.  Note that either upper or lower case may be
!            used.  A more detailed explanation is provided in the
!            portion of the code where these strings are declared.
! ------------------------------
!  The following three strings designate the type of sounding data.
!
!  if height(1:1)='h' or 'H', the sounding is given on height levels
!  if height(1:1)='p' or 'P', the sounding is given on pressure levels
!
!  if therm(1:1) ='t' or 'T', the sounding is specified in temperature
!  if therm(1:1) ='p' or 'P', the sounding is specified in potential
!                             temperature.
!
!  if humid(1:1) ='s' or 'S', the soundings uses specific humidity
!  if humid(1:1) ='r' or 'R', the sounding uses relative humidity,
!  if humid(1:1) ='d' or 'D', the soundings uses dewpoint temperature,
!
!  if wind(1:1) ='x' or 'Y', the sounding is specified in x and y
!                            component of velocity (u and v)
!  if wind(1:1) ='d' or 'D', the sounding is specified in direction
!                            and speed in m/s.
!  if wind(1:1) ='k' or 'K', the sounding is specified in direction
!                            and speed in knots.
!--------------------------------
!  Record 7: Ground-level height (m) and the correspoding pressure
!            (Pascal) when sounding is specified at height levels.
!            When it is given on pressure levels, this record is
!            not used. The last sounding data is assumed to be at
!            the ground level (z=0).
!  Record 8: number of data levels in sounding (variable lvlsnd )
!  Record 9: a line of data column labels (not read in)
!  Record 10 to the end:
!            sounding data in the order 'z/p, pt/t, qv/rh, u, v'
!
!  For records 10  to the end, there is one line of data
!            corresponding to each sounding level.
!
!  Important convention:
!            Line 10 corresponds to the level k = lvlsnd
!            Line 11 corresponds to the level k = (lvlsnd -1)
!            etc.
!
!  Units of the data:
!            pressure: Pascals, height: meters, temperature, dewpoint,
!            and potential temperature: degrees Kelvin, specific
!            humidity kg/kg, relative humidity: 0 to 1.
!
!  CAUTION: The sounding data have to be listed in the order of
!           decreasing height or increasing pressure.
!
!-----------------------------------------------------------------------
!
      WRITE(*,*)  year,month,day,hour,minute,second
      WRITE(sndtime,'(i2,a,i2,a,a)') ihr,':',imin,':','00 CST'
      IF(ihr <=9 ) sndtime(1:1)='0'
      IF(imin <=9 ) sndtime(4:4)='0'
      WRITE(snddate,'(i3,i3,a,i5,a,f6.1,a,f6.1,a)')  iday,imo,',',year, & 
        ' - Domain speed plused (umove=',ustorm,',vmove=',ustorm,')'
      WRITE(sndlocation,'(a,f8.2,3x,a,f8.2,3x,a,f7.1)')                 &
       'SLAT=',slat(isnd),'SLON=',slon(isnd),'SELV=',selev(isnd)

      OPEN(3,FILE=trim(ofile(isnd))//'forARPS',STATUS='unknown')
      WRITE(3,*) '1-D Sounding Input for ARPS'
      WRITE(3,*) 'supercell storm'
      WRITE(3,'(1x,a)')  sndtime
      WRITE(3,'(a)')  snddate
      WRITE(3,'(1x,a)')  sndlocation
      WRITE(3,*) '''height'' ''potential temperature'' ''specific humidity'' ''uv'' '
      WRITE(3,*)   shght(1,isnd), spres(1,isnd)*100
      WRITE(3,*)   nlevs
      WRITE(3,*)   ' height potential  temperature   SH    U         V' 

!
      DO k=nlevs,1,-1
        WRITE(3,'(1x,6(F14.5))')                                        &
              shght(k,isnd),                                            &
              stheta(k,isnd),sqv(k,isnd),                               &
              su(k,isnd),sv(k,isnd)
      END DO
      CLOSE(3)
      write(6,'(a,a,a)')                                                &
      'Sounding file ',trim(ofile(isnd))//'forARPS',' was produced.'

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

    END DO  ! the number of sounding
  ELSE
    WRITE(6,'(a)') ' Error reading data.  EXTSND ends'
  END IF
  STOP

END PROGRAM ARPSEXTSND 
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE COLINT                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE colint(nx,ny,nz,maxsnd,nsnd,                                 & 5,8
           xs,ys,zp,xpt,ypt,ipt,jpt,                                    &
           uprt, vprt, ptprt, pprt, qvprt,                              &
           ubar, vbar, ptbar, pbar, qvbar,                              &
           su,sv,stheta,spres,shght,sqv,selev,                          &
           dxfld,dyfld,rdxfld,rdyfld,                                   &
           slopey,alphay,betay,                                         &
           tem1,nlevs)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolates ARPS history data in the horizontal to create
!  a column of data located at point xpt, ypt.
!
!  Bilinear interpolation is used.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  April 1992.
!
!  MODIFICATION HISTORY:
!
!  October, 1992 (K. Brewster)
!  Conversion to ARPS 3.0.
!
!  October, 1994 (K. Brewster)
!  Conversion to ARPS 4.0.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx,ny,nz Dimensions of ARPS grids.
!
!    xs       x coordinate of scalar points in physical/comp. space (m)
!    ys       y coordinate of scalar points in physical/comp. space (m)
!    zp       z coordinate of scalar grid points in physical space (m)
!
!    xpt      x coordinate of desired sounding (m)
!    ypt      y coordinate of desired sounding (m)
!
!    ipt      i index of grid point just west of xpt,ypt
!    jpt      j index of grid point just south of xpt,ypt
!
!    uprt     x component of perturbation velocity (m/s)
!    vprt     y component of perturbation velocity (m/s)
!
!    ptprt    Perturbation potential temperature (K)
!    pprt     Perturbation pressure (Pascal)
!
!    qvprt    Perturbation water vapor mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!  OUTPUT:
!
!    su       Interpolated u wind component.  (m/s)
!    sv       Interpolated v wind component.  (m/s)
!    stheta   Interpolated potential temperature (K).
!    spres    Interpolated pressure. (Pascals)
!    shght    Interpolated height (meters)
!    sqv      Interpolated water vapor mixing ratio (kg/kg).
!    selev    Interpolated surface elevation (m)
!    nlevs    Number of above-ground sounding levels.
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
!  Arguments -- location data
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Dimensions of ARPS grids.
  INTEGER :: maxsnd
  INTEGER :: nsnd
  REAL :: xs(nx)               ! x coordinate of grid points in
                               ! physical/comp. space (m)
  REAL :: ys(ny)               ! y coordinate of grid points in
                               ! physical/comp. space (m)
  REAL :: zp(nx,ny,nz)         ! z coordinate of grid points in
                               ! physical space (m)
  REAL :: xpt(maxsnd)          ! location to find in x coordinate (m)
  REAL :: ypt(maxsnd)          ! location to find in y coordinate (m)
  INTEGER :: ipt(maxsnd)       ! i index to the west of desired
                               ! location
  INTEGER :: jpt(maxsnd)       ! j index to the south of desired
                               ! location
!
!-----------------------------------------------------------------------
!
!  Arguments -- model data
!
!-----------------------------------------------------------------------
!
  REAL :: uprt   (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: vprt   (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: ptprt  (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: qvprt  (nx,ny,nz)    ! Perturbation water vapor specific
                               ! humidity (kg/kg)

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

  REAL :: dxfld(nx)
  REAL :: dyfld(ny)
  REAL :: rdxfld(nx)
  REAL :: rdyfld(ny)
  REAL :: slopey(nx,ny,nz)
  REAL :: alphay(nx,ny,nz)
  REAL :: betay(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Arguments -- Extracted sounding variables
!
!-----------------------------------------------------------------------
!
  REAL :: su(nz,maxsnd)
  REAL :: sv(nz,maxsnd)
  REAL :: stheta(nz,maxsnd)
  REAL :: sqv(nz,maxsnd)
  REAL :: spres(nz,maxsnd)
  REAL :: shght(nz,maxsnd)
  REAL :: selev(maxsnd)
  INTEGER :: nlevs
!
!-----------------------------------------------------------------------
!
!  Temporary work arrays
!
!-----------------------------------------------------------------------
!
  REAL :: tem1(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Functions called
!
!-----------------------------------------------------------------------
!
  REAL :: aint2d
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,kk,isnd
  REAL :: w2,w3
  REAL :: t1,t2,t3,hmid,tmid,qvsat,rh
  INTEGER :: iorder
  PARAMETER (iorder = 2)
  REAL :: pntint2d
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$*  inline routine (f_qvsat)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nlevs=nz-2

  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,zp,                                         &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      shght(k,isnd)=pntint2d(nx,ny,                                     &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),zp(1,1,k),                           &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k)=ptbar(i,j,k)+ptprt(i,j,k)
      END DO
    END DO
  END DO
  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,tem1,                                       &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      stheta(k,isnd)=pntint2d(nx,ny,                                    &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),tem1(1,1,k),                         &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k)=pbar(i,j,k)+pprt(i,j,k)
      END DO
    END DO
  END DO
  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,tem1,                                       &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      spres(k,isnd)=pntint2d(nx,ny,                                     &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),tem1(1,1,k),                         &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k)=qvbar(i,j,k)+qvprt(i,j,k)
      END DO
    END DO
  END DO
  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,tem1,                                       &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      sqv(k,isnd)=pntint2d(nx,ny,                                       &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),tem1(1,1,k),                         &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Get height at scalar points, since zp was defined at w points.
!
!-----------------------------------------------------------------------
!
  DO isnd=1,nsnd
    selev(isnd)=shght(2,isnd)
    DO k=1,nz-1
      shght(k,isnd)=0.5*(shght(k+1,isnd)+shght(k,isnd))
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Form total u wind component at scalar points
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k)=0.5*(ubar(i,j,k)+ubar(i+1,j,k))+                    &
                    0.5*(uprt(i,j,k)+uprt(i+1,j,k))
      END DO
    END DO
  END DO

  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,tem1,                                       &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      su(k,isnd)=pntint2d(nx,ny,                                        &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),tem1(1,1,k),                         &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Form total v wind component at scalar points
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k)=0.5*(vbar(i,j,k)+vbar(i,j+1,k)) +                   &
                    0.5*(vprt(i,j,k)+vprt(i,j+1,k))
      END DO
    END DO
  END DO

  CALL setdrvy(nx,ny,nz,                                                &
               1,nx-1,1,ny-1,1,nz-1,                                    &
               dyfld,rdyfld,tem1,                                       &
               slopey,alphay,betay)
  DO isnd=1,nsnd
    DO k=2,nz-1
      sv(k,isnd)=pntint2d(nx,ny,                                        &
               1,nx-1,1,ny-1,                                           &
               iorder,xs,ys,xpt(isnd),ypt(isnd),                        &
               ipt(isnd),jpt(isnd),tem1(1,1,k),                         &
               dxfld,dyfld,rdxfld,rdyfld,                               &
               slopey(1,1,k),alphay(1,1,k),betay(1,1,k))
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Get a value at the surface, by extrapolating from the
!  2nd and third levels.
!
!-----------------------------------------------------------------------
!
  DO isnd=1,nsnd
    shght(1,isnd)=selev(isnd)
    w3=(shght(1,isnd)-shght(2,isnd))                                    &
        /(shght(3,isnd)-shght(2,isnd))
    w2=(1.-w3)
    su(1,isnd)=w2*    su(2,isnd) + w3*    su(3,isnd)
    sv(1,isnd)=w2*    sv(2,isnd) + w3*    sv(3,isnd)
    IF(stheta(3,isnd) > stheta(2,isnd)) THEN
      stheta(1,isnd)=w2*stheta(2,isnd) + w3*stheta(3,isnd)
    ELSE
      stheta(1,isnd)=stheta(2,isnd)
    END IF
!
!-----------------------------------------------------------------------
!
!  Integrate downward to get the pressure at level 1.
!
!-----------------------------------------------------------------------
!
    t3=stheta(3,isnd)*(spres(3,isnd)/100000.)**rddcp
    t2=stheta(2,isnd)*(spres(2,isnd)/100000.)**rddcp
    hmid=0.5*(shght(2,isnd)+shght(1,isnd))
    tmid=t3+((shght(3,isnd)-hmid)/                                      &
            (shght(3,isnd)-shght(2,isnd)))*(t2-t3)
    spres(1,isnd)=spres(2,isnd)*                                        &
           EXP(g*(shght(2,isnd)-shght(1,isnd))/(rd*tmid))
!
!-----------------------------------------------------------------------
!
!  Use constant RH to extrapolate qv to level 1.
!
!-----------------------------------------------------------------------
!
    qvsat = f_qvsat( spres(2,isnd), t2 )     ! saturated S.H.
    rh=AMIN1((sqv(2,isnd)/qvsat),1.)         ! R.H.
    t1=stheta(1,isnd)*(spres(1,isnd)/100000.)**rddcp
    qvsat = f_qvsat( spres(1,isnd), t1 )
    sqv(1,isnd)=rh*qvsat
  END DO
  RETURN
END SUBROUTINE colint
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE CNVSND                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

SUBROUTINE cnvsnd(su,sv,stheta,spres,sqv,slon,                          & 1,1
           sdrct,ssped,stemp,sdewp,nlevs)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Converts units of data extracted from ARPS history data to
!  those required of sounding data. Determines direction and
!  speed from u and v wind components.
!
!  Dew-point formula from Bolton, 1980, MWR pp 1046-1053.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  April 1992.
!
!  MODIFICATION HISTORY:
!
!  October, 1992 (K. Brewster)
!  Conversion to ARPS 3.0.
!
!  10/28/1992 (K. Brewster)
!  Special allowance for low qv-to-dew pt
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    su       Sounding u wind component.  (m/s)
!    sv       Sounding v wind component.  (m/s)
!    stheta   Sounding potential temperature (K).
!    spres    Sounding pressure. (Pascals)
!    sqv      Sounding specific humidity
!    slon     Sounding longitude
!    nlevs    Number of above-ground sounding levels.
!
!  OUTPUT:
!
!    spres    Sounding pressure. (mb)
!    sdrct    Sounding wind direction (degrees from north)
!    ssped    Sounding wind speed (m/s)
!    stemp    Sounding temperature (degrees C)
!    sdewp    Sounding dew point temperature (degrees C)
!
!-----------------------------------------------------------------------
!
!  Variable declarations
!
!-----------------------------------------------------------------------
!
!  Input arguments
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nlevs             ! Number of above-ground sounding levels
  REAL :: su    (nlevs)        ! Sounding u wind component (m/s)
  REAL :: sv    (nlevs)        ! Sounding v wind component (m/s)
  REAL :: stheta(nlevs)        ! Sounding potential temperature (K)
  REAL :: spres (nlevs)        ! Sounding pressure. (Pascals)
  REAL :: sqv   (nlevs)        ! Sounding specific humidity (g/g)
  REAL :: slon                 ! Sounding longitude (degrees E)
!
!-----------------------------------------------------------------------
!
!  Output arguments
!
!-----------------------------------------------------------------------
!
  REAL :: sdrct(nlevs)         ! Sounding wind direction
                               ! (degrees from north)
  REAL :: ssped(nlevs)         ! Sounding wind speed (m/s)
  REAL :: stemp(nlevs)         ! Sounding temperature (degrees C)
  REAL :: sdewp(nlevs)         ! Sounding dew point temperature
                               ! (degrees C)
!
!-----------------------------------------------------------------------
!
!  Constants
!
!-----------------------------------------------------------------------
!
  REAL :: rd2dg
  PARAMETER (rd2dg=(180./3.141592654))
  REAL :: rhmax
  PARAMETER (rhmax=1.0)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: k
  REAL :: smix,e,bige,alge
!
  DO k=1,nlevs
!
!-----------------------------------------------------------------------
!
!  Convert u,v to direction and speed
!
!-----------------------------------------------------------------------
!
    CALL uvrotdd(1,1,slon,su(k),sv(k),sdrct(k),ssped(k))
!
!-----------------------------------------------------------------------
!
!  Convert pressure from Pascals to mb
!
!-----------------------------------------------------------------------
!
    spres(k)=spres(k)*0.01
!
!-----------------------------------------------------------------------
!
!  Convert theta to temperature in degrees C.
!
!-----------------------------------------------------------------------
!
    stemp(k)=stheta(k)*((spres(k)/1000.)**rddcp)
    stemp(k)=stemp(k)-273.15
!
!-----------------------------------------------------------------------
!
!  Convert specific humidity to dew-point in degrees C.
!
!-----------------------------------------------------------------------
!
    IF( sqv(k) > 0.) THEN
      smix=sqv(k)/(1.-sqv(k))
      e=(spres(k)*smix)/(0.62197 + smix)
      bige=e/( 1.001 + ( (spres(k) - 100.) / 900.) * 0.0034)
      alge = ALOG(bige/6.112)
      sdewp(k)= (alge * 243.5) / (17.67 - alge)
    ELSE
      sdewp(k)= stemp(k)-30.
    END IF
  END DO
!
  RETURN 

END SUBROUTINE cnvsnd