!
!##################################################################
!##################################################################
!######                                                      ######
!######                   PROGRAM EXT2ARPS                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


PROGRAM ext2arps,172
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Converts external forecast files coordinates and variables
!  to those used in ARPS and writes the information in a
!  history dump file.
!
!  The actual reading of the external files is done by
!  subroutine rdextfil.  In many cases, it will only be necessary
!  to make changes in extdims.inc and rdextfil to customize the
!  processing for the user's specific data source.
!
!  This program could become more general in the future, but
!  for now we assume incoming data are supplied (or converted to)
!  pressure(Pa), height(m), temperature(K), specific humidty (kg/kg),
!  and u,v (m/s).  Units should be converted from their original
!  forms to these in rdextfil.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  Sept, 1994  after MAPS2ARPS
!
!  MODIFICATION HISTORY:
!  Oct, 1994 (KB)
!  Changed input times to be more flexible.
!
!  17 Jan, 1995 (KB)
!  Modified processing so that mean sounding would be for constant
!  pressure surfaces rather than constant ARPS levels.
!
!  23 Jan, 1995 (KB)
!  Added temporary array avgzs_ext, replacing the use of tem1 which
!  would not work if number of ext levels is greater than ARPS.
!
!  9 August, 1995 (KB)
!  Modified creation of mean soundings and storage of external
!  pressure arrays to allow for incoming data on other than
!  pressure vertical coordinates.  Added switch for buoyancy
!  balancing of vertical pressure gradient.
!
!  26 April, 1996 (KB)
!  Version 2.0, including:  Corrects divergence error
!  using a function either linear in k-level, physical height or
!  potential temperature.  Also, moisture is interpolated and
!  extrapolated using RHstar [ RHstar=sqrt(1. - RH) ], which should
!  be more linear in height, rather than qv.   Extrapolation above
!  and below available model data are constant-in-height perturbation
!  fields for all variables except assigned a standard lapse rate near
!  the ground and dT/dz=0 at the top.
!
!  30 April 1996 (KB)
!  Final preparation of version 2.0, added O'Brien option and
!  order of interpolation option to input namelist.
!
!  8 November 1996 (KB)
!  Version 2.1
!  Changed the polynomial interpolation scheme for better
!  accuracy and eliminated chance for a discontinuity at external
!  grid cell boundaries.  Linear or quadratic interpolation
!  using local polynomial fit are supported (iorder=1 or 2).
!
!  3 March 1996 (KB)
!  Pressure is now interpolated with polynomials of pressure
!  rather than ln(p).
!
!  10 June 1997 (Fanyou Kong -- CMRP)
!  Include surface 2D data analysis from external data sets, a new
!  subroutine MKARPS2D is then added in file 'extlib23.f'.
!
!  16 July 1997 (Zonghui Huo)
!  Added 'getcoamps.f' to read the COAMPS data on sigma-z levels.
!
!  06/16/1998 (Dennis Moon, SSESCO)
!  Added support for RUC 2 on the AWIPS 211 grid. extdopt=7
!
!  06/16/1998 (Dennis Moon, SSESCO)
!  Added support for NCEP global re-analysis on T62 Gaussian lat/lon
!  grid. This does not work for ARPS grids which cross the 0 degree
!  meridian. extdopt=8
!
!  29 Oct 1998 (R. Carpenter, G. Bassett)
!  Added RUC-2 GEMPAK & ETA GEMPAK grid #104.  Four digit years and
!  extdname now used in constructing filenames of GEMPAK files.
!
!  15 Sep 1999 (K. Brewster)
!  Added the microphysical variables to the dtadump calls so they
!  will be used in those cases when available.  Added check for
!  positive definiteness to microphysical variables.
!
!  2000/03/23 (Donghai Wang)
!  Added NCEP AVN GRIB global data, grid #3.
!
!  2000/07/28 (Ming Xue)
!  Converted to F90 free format and implemented dynamic memory
!  allocation.
!
!  2000/08/14 (Gene Bassett)
!  Added multiple soil types, sfcdata variables and grid staggering
!  for use with arps history format of external model data.  Added options
!  allowing ext2arps to produce results similar to arpsintrp (intropt,
!  ext_lbc, ext_vbc).
!
!  2001/05/31 (Gene Bassett)
!  Added exttrnopt, an option to use input grid terrain for output grid.
!  Also corrected a bug introduced when updating the use of arps as an 
!  external model (2000/08) which caused only the first external file 
!  to have u&v rotated to arps map projection.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Input and output grid dimensions
!
!-----------------------------------------------------------------------

  INTEGER :: nx ! Number of grid points in the x-dir of output arps grid
  INTEGER :: ny ! Number of grid points in the y-dir of output arps grid
  INTEGER :: nz ! Number of grid points in the z-dir of output arps grid
!
  INTEGER :: nx_ext, ny_ext, nz_ext ! dimensions of external data grid

  INTEGER :: nstyps
!
!-----------------------------------------------------------------------
!
!  Space for mean sounding
!
!-----------------------------------------------------------------------
!
  INTEGER :: lvlprof
  REAL :: depthp

  PARAMETER( lvlprof=601, depthp = 3.0E4 )
!
!-----------------------------------------------------------------------
!
!  lvlprof:
!
!  The ARPS interpolates unevenly-spaced data from the external
!  data set to evenly-spaced data for use in defining the base
!  state atmosphere.  In this process, an intermediate sounding is
!  generated at evenly-spaced altitudes, with the accuracy of the
!  associated interpolation controlled by the parameter lvlprof.
!  The larger lvlprof, the more accurate the interpolation
!  (we recommend using lvlprof=200 for a model run with about 50
!  points in the vertical direction).  Using the intermediate
!  sounding, the ARPS then generates a base state model sounding
!  for the particular vertical grid resolution chosen (i.e.,
!  the number of points in the vertical, nz, and the vertical grid
!  spacing, dz).
!
!  depthp:
!
!  The depth of atmosphere over which the interpolated profiles
!  will be defined.  depthp should be greater than or equal to
!  (nz-3)*dz, i.e., larger than the physical depth of the model
!  domain.  Otherwise, the code will extrapolate for gridpoints
!  outside the domain, leading to possible inconsistencies.  At all
!  costs, any such extrapolation should be avoided.
!
!-----------------------------------------------------------------------
!
  REAL :: psnd(lvlprof),                                                &
       zsnd(lvlprof),                                                   &
       tsnd(lvlprof),                                                   &
       ptsnd(lvlprof),                                                  &
       rhssnd(lvlprof),                                                 &
       qvsnd(lvlprof),                                                  &
       rhosnd(lvlprof),                                                 &
       usnd(lvlprof),                                                   &
       vsnd(lvlprof),                                                   &
       dumsnd(lvlprof),                                                 &
       plsnd(lvlprof)
!
!-----------------------------------------------------------------------
!
!  ARPS grid variables
!
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Vertical component of Cartesian velocity at a given
!             time level (m/s)
!    ptprt    Perturbation potential temperature at a given time
!             level (K)
!    pprt     Perturbation pressure at  a given time level (Pascal)
!    qv       Water vapor specific humidity at a given time level (kg/kg)
!    qc       Cloud water mixing ratio at a given time level (kg/kg)
!    qr       Rainwater mixing ratio at a given time level (kg/kg)
!    qi       Cloud ice mixing ratio at a given time level (kg/kg)
!    qs       Snow mixing ratio at a given time level (kg/kg)
!    qh       Hail mixing ratio at a given time level (kg/kg)
!
!    ubar     Base state x-velocity component (m/s)
!    vbar     Base state y-velocity component (m/s)
!    wbar     Base state vertical velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    hterain  Terrain height (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    tsfc     Temperature at ground (K) (in top 1 cm layer)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture in the top 1 cm layer
!    wetdp    Deep soil moisture in the deep 1 m layer
!    wetcanp  Canopy water amount
!
!-----------------------------------------------------------------------
!
  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 :: xscl(:)    ! The x-coordinate of scalar points.
  REAL, allocatable :: yscl(:)    ! The y-coordinate of scalar points.
  REAL, allocatable :: zp(:,:,:)  ! The physical height coordinate defined at
                                  ! w-point of the staggered grid.
  REAL, allocatable :: zs(:,:,:)  ! The physical height of scalar points.
  REAL, allocatable :: j1(:,:,:)  ! Coordinate transformation Jacobian -d(zp)/dx.
  REAL, allocatable :: j2(:,:,:)  ! Coordinate transformation Jacobian -d(zp)/dy.
  REAL, allocatable :: j3(:,:,:)  ! Coordinate transformation Jacobian  d(zp)/dz.
  REAL, allocatable :: aj3z(:,:,:)! Coordinate transformation Jacobian defined
                                  ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL, allocatable :: hterain(:,:) ! The height of the terrain. (m)
!wdt update
  REAL, allocatable :: htrn_ext(:,:) ! The height of the terrain. (m)
  REAL, allocatable :: mapfct(:,:,:)! Map factor at scalar, u, and v points

  REAL, allocatable :: u(:,:,:)     ! Total u-velocity (m/s)
  REAL, allocatable :: v(:,:,:)     ! Total v-velocity (m/s)
  REAL, allocatable :: w(:,:,:)     ! Total w-velocity (m/s)
  REAL, allocatable :: pprt(:,:,:)  ! Perturbation pressure from that
                                    ! of base state atmosphere (Pascal)
  REAL, allocatable :: ptprt(:,:,:) ! Perturbation potential temperature
                                    ! from that of base state atmosphere (K)
!
  REAL, allocatable :: qv(:,:,:)    ! 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 :: pbar(:,:,:)  ! Base state pressure (Pascal)
  REAL, allocatable :: ptbar(:,:,:) ! Base state potential temperature (K)
  REAL, allocatable :: qvbar(:,:,:) ! Base state water vapor specific humidity
                                    ! (kg/kg)
  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 :: rhobar(:,:,:)! Base state density (kg/m3).
  REAL, allocatable :: rhostr(:,:,:)! Base state density rhobar times j3.
  REAL, allocatable :: wcont(:,:,:) ! Contravariant vertical velocity (m/s)

  REAL, allocatable :: csndsq(:,:,:)! Speed of sound squared (m**2/s**2)


  REAL, allocatable :: tsfc(:,:,:)   ! Ground surface temperature (K)
  REAL, allocatable :: tsoil(:,:,:)  ! Deep soil temperature (K)
  REAL, allocatable :: wetsfc(:,:,:) ! Surface soil moisture
  REAL, allocatable :: wetdp(:,:,:)  ! Deep soil moisture
  REAL, allocatable :: wetcanp(:,:,:)! Water amount on canopy
  REAL, allocatable :: snowdpth(:,:) ! Snow depth (m)

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

  REAL, allocatable :: tem1(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem2(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem3(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem4(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem5(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem6(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem7(:,:,:)   ! Temporary work array.
  REAL, allocatable :: tem8(:,:,:)   ! Temporary work array.

  REAL, allocatable :: xs2d(:,:)
  REAL, allocatable :: ys2d(:,:)
  REAL, allocatable :: xu2d(:,:)
  REAL, allocatable :: yu2d(:,:)
  REAL, allocatable :: xv2d(:,:)
  REAL, allocatable :: yv2d(:,:)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: xw & yw
  REAL, allocatable :: xw(:,:)
  REAL, allocatable :: yw(:,:)

  INTEGER, allocatable :: iscl(:,:) ! i index of scalar point
                                    ! in external array.
  INTEGER, allocatable :: jscl(:,:) ! j index of scalar point
                                    ! in external array.
  INTEGER, allocatable :: iu(:,:)   ! i index of u-velocity point
                                    ! in external array.
  INTEGER, allocatable :: ju(:,:)   ! j index of u-velocity point
                                    ! in external array.
  INTEGER, allocatable :: iv(:,:)   ! i index of v-velocity point
                                    ! in external array.
  INTEGER, allocatable :: jv(:,:)   ! j index of v-velocity point
                                    ! in external array.

  INTEGER :: iproj_arps           ! ARPS map projection (tem storage)
  REAL :: scale_arps              ! ARPS map scale (tem storage)
  REAL :: trlon_arps              ! ARPS true-longitude (tem storage)
  REAL :: latnot_arps(2)          ! ARPS true-latitude(s) (tem storage)
  REAL :: xorig_arps,yorig_arps   ! ARPS map projection origin (tem storage)
!
!-----------------------------------------------------------------------
!
!  NAMELIST parameter (In arps.input)
!
!-----------------------------------------------------------------------
!
  INCLUDE 'ext2arps.inc'
!
!-----------------------------------------------------------------------
!
!  External forecast variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext            ! external data map projection
  REAL :: scale_ext               ! external data map scale factor
  REAL :: trlon_ext               ! external data true longitude
  REAL :: latnot_ext(2)           ! external data true latitude(s)
  REAL :: x0_ext,y0_ext           ! external data origin

  REAL, allocatable :: x_ext(:)       ! external data x-coordinate
  REAL, allocatable :: y_ext(:)       ! external data y-coordinate
  REAL, allocatable :: xu_ext(:)       ! external data u x-coordinate
  REAL, allocatable :: yu_ext(:)       ! external data u y-coordinate
  REAL, allocatable :: xv_ext(:)       ! external data v x-coordinate
  REAL, allocatable :: yv_ext(:)       ! external data v y-coordinate
  REAL, allocatable :: lat_ext(:,:)    ! external data latidude
  REAL, allocatable :: lon_ext(:,:)    ! external data longitude
  REAL, allocatable :: latu_ext(:,:)   ! external data latidude (x-stag)
  REAL, allocatable :: lonu_ext(:,:)   ! external data longitude (x-stag)
  REAL, allocatable :: latv_ext(:,:)   ! external data latidude (y-stag)
  REAL, allocatable :: lonv_ext(:,:)   ! external data longitude (y-stag)
  REAL, allocatable :: zs_ext(:,:,:)   ! external data physical height (m)

  REAL, allocatable :: p_ext(:,:,:)    ! Pressure (Pascals)
  REAL, allocatable :: hgt_ext(:,:,:)  ! Height (m)
  REAL, allocatable :: zp2_ext(:,:,:)  ! external w physical height (m)
  REAL, allocatable :: zp_ext(:,:,:)   ! Height (m) (on arps grid)
  REAL, allocatable :: t_ext(:,:,:)    ! Temperature (K)
  REAL, allocatable :: u_ext(:,:,:)    ! Eastward wind component
  REAL, allocatable :: v_ext(:,:,:)    ! Northward wind component
  REAL, allocatable :: w_ext(:,:,:)    ! Vertical wind component
  REAL, allocatable :: qv_ext(:,:,:)   ! Specific humidity (kg/kg)
  REAL, allocatable :: qc_ext(:,:,:)   ! Cloud H2O mixing ratio (kg/kg)
  REAL, allocatable :: qr_ext(:,:,:)   ! Rain  H2O mixing ratio (kg/kg)
  REAL, allocatable :: qi_ext(:,:,:)   ! Ice   H2O mixing ratio (kg/kg)
  REAL, allocatable :: qs_ext(:,:,:)   ! Snow  H2O mixing ratio (kg/kg)
  REAL, allocatable :: qh_ext(:,:,:)   ! Hail  H2O mixing ratio (kg/kg)

  REAL, allocatable :: tsfc_ext   (:,:,:)   ! Temperature at surface (K)
  REAL, allocatable :: tsoil_ext  (:,:,:)   ! Deep soil temperature (K)
  REAL, allocatable :: wetsfc_ext (:,:,:)   ! Surface soil moisture
  REAL, allocatable :: wetdp_ext  (:,:,:)   ! Deep soil moisture
  REAL, allocatable :: wetcanp_ext(:,:,:)   ! Canopy water amount
  REAL, allocatable :: snowdpth_ext(:,:)  ! Snow depth (m)

  INTEGER, allocatable :: soiltyp_ext (:,:,:)! Soil type
  REAL, allocatable ::    stypfrct_ext(:,:,:)! Soil type fraction
  INTEGER, allocatable :: vegtyp_ext (:,:)   ! Vegetation type
  REAL, allocatable ::    lai_ext    (:,:)   ! Leaf Area Index
  REAL, allocatable ::    roufns_ext (:,:)   ! Surface roughness
  REAL, allocatable ::    veg_ext    (:,:)   ! Vegetation fraction

  REAL, allocatable :: trn_ext    (:,:)   ! External terrain (m)
  REAL, allocatable :: psfc_ext   (:,:)   ! Surface pressure (Pa)

  REAL, allocatable :: t_2m_ext (:,:)
  REAL, allocatable :: rh_2m_ext(:,:)
  REAL, allocatable :: u_10m_ext(:,:)
  REAL, allocatable :: v_10m_ext(:,:)

  REAL, allocatable :: dxfld(:)
  REAL, allocatable :: dyfld(:)
  REAL, allocatable :: rdxfld(:)
  REAL, allocatable :: rdyfld(:)
  REAL, allocatable :: dxfldu(:)
  REAL, allocatable :: dyfldu(:)
  REAL, allocatable :: rdxfldu(:)
  REAL, allocatable :: rdyfldu(:)
  REAL, allocatable :: dxfldv(:)
  REAL, allocatable :: dyfldv(:)
  REAL, allocatable :: rdxfldv(:)
  REAL, allocatable :: rdyfldv(:)
  REAL, allocatable :: tem1_ext(:,:,:)    ! Temporary work array
  REAL, allocatable :: tem2_ext(:,:,:)    ! Temporary work array
  REAL, allocatable :: tem3_ext(:,:,:)    ! Temporary work array
  REAL, allocatable :: tem4_ext(:,:,:)    ! Temporary work array
  REAL, allocatable :: tem5_ext(:,:,:)    ! Temporary work array

  REAL, allocatable :: xa_ext(:,:)
  REAL, allocatable :: ya_ext(:,:)
  REAL, allocatable :: avgzs_ext(:,:,:)
!
!-----------------------------------------------------------------------
!
!  include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'adjust.inc'
!
!-----------------------------------------------------------------------
!
!  Non-dimensional smoothing coefficient
!
!-----------------------------------------------------------------------
!
  REAL :: smfct1,smfct2
  PARAMETER( smfct1=0.5, smfct2=-0.5 )
!
  REAL :: rhmax
  PARAMETER ( rhmax = 1.0)
!
!-----------------------------------------------------------------------
!
!  Latitude and longitude for some diagnostic printing,
!  e.g. to compare to an observed sounding
!
!-----------------------------------------------------------------------
!
  REAL :: latdiag,londiag
  PARAMETER (latdiag=34.5606,londiag=-103.0820)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: basdmpfn
  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfinit
  CHARACTER (LEN=9) :: julfname
!
  INTEGER :: i,j,k,ksmth,istatus
  INTEGER :: iyr,imo,iday,ihr,imin,isec,jldy
  INTEGER :: ifhr,ifmin,ifsec,mfhr
  INTEGER :: myr,initsec,iabssec,jabssec,kftime
  INTEGER :: ifile,iprtopt,iniotfu,lbasdmpf,onvf,grdbas
  INTEGER :: iextmn,iextmx,jextmn,jextmx
  INTEGER :: idiag,jdiag

  REAL :: latnot(2)
  REAL :: amin,amax
  REAL :: qvmin,qvmax,qvval
  REAL :: csconst,pconst
  REAL :: deltaz,tv_ext
  REAL :: pres,temp,qvsat,rh,tvbar,qvprt,qtot
  REAL :: xdiag,ydiag,dd,dmin,latd,lond
  REAL :: ppasc,pmb,tc,tdc,theta,smix,e,bige,alge,dir,spd

!wdt update
  CHARACTER (LEN=80) :: soiloutfl,temchar,ternfn
  CHARACTER (LEN=80) :: timsnd
  INTEGER :: lfn, tmstrln, lternfn
  INTEGER :: tsfcout,tsoilout,wetsout,wetdout,wetcout,snowdout
  INTEGER :: isnow,jsnow,ii,jj
  REAL :: xumin,xumax,yvmin,yvmax

  INTEGER :: strlen, ireturn
  CHARACTER (LEN=3) :: FMT
  CHARACTER (LEN=100) :: tmp_ch

  INTEGER :: is,nstyp_ext

!wdt-williams 2002-01-17 GMB ntmerge
  DOUBLE PRECISION :: ntmergeinv, mfac
  INTEGER :: idist
!
!-----------------------------------------------------------------------
!
!  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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  mgrid = 1
  nestgrd = 0

  nstyps = 4
  nstyp_ext = 0

  DO k=1,lvlprof
    dumsnd(k) = 0.0
  END DO
!
!-----------------------------------------------------------------------
!
!  Call initpara to read in ARPS parameters of namelists
!
!-----------------------------------------------------------------------
!
  CALL initpara(nx,ny,nz,nstyps)

  CALL initadas ! The adjust namelist block is used by ext2arps.  Other
                ! adas parameters are not used by this program,
                ! but the name list blocks need to be read in sequence on 
                ! Cray machines such as J90. 

  myr=MOD(year,100)
  WRITE(julfinit,'(i2.2,i3.3,i2.2,i2.2)') myr,jday,hour,minute
  CALL ctim2abss(year,month,day,hour,minute,second,initsec)
!
!-----------------------------------------------------------------------
!
!  Read in additional namelists for external file specifications.
!
!-----------------------------------------------------------------------
!
  extdopt  = 0
  extdfmt  = 10
  dir_extd = './'
  extdname = 'may20'
  nextdfil = 1
  extdtime(1) = '1977-05-20.21:00:00+000:00:00'
  iorder   = 2
  intropt  = 1
  nsmooth  = 1
  ext_lbc  = 1
  ext_vbc  = 1
  exttrnopt = 0
!wdt-williams 2002-01-17 GMB ntmerge
  extntmrg = 1

  READ (5,extdfile)

  WRITE (6, '(/1x,a)')    '&extdfile'

  strlen = LEN( dir_extd )
  CALL strlnth( dir_extd, strlen)
  WRITE (6, '(3x,a,a,a)')                                               &
                       'dir_extd   = ''', dir_extd(1:strlen), ''','

  strlen = LEN( extdname )
  CALL strlnth( extdname, strlen)
  WRITE (6, '(3x,a,a,a)')                                               &
                       'extdname   = ''', extdname(1:strlen), ''','

  WRITE (6, '(3x,a,i4,a)')    'extdopt   = ', extdopt,  ','
  WRITE (6, '(3x,a,i4,a)')    'extdfmt   = ', extdfmt,  ','
  WRITE (6, '(3x,a,i4,a)')    'nextdfil  = ', nextdfil, ','

  DO i=1,nextdfil
    strlen = LEN( extdtime(i) )
    CALL strlnth( extdtime(i), strlen)
    WRITE (6, '(3x,a,i2.2,a,a,a)')                                      &
                       'extdtime(',i,') = ''', extdtime(i), ''','
  END DO

  WRITE (6, '(3x,a,i4,a)')    'iorder    = ', iorder,   ','
  WRITE (6, '(3x,a,i4,a)')    'intropt   = ', intropt,  ','
  WRITE (6, '(3x,a,i4,a)')    'nsmooth   = ', nsmooth,  ','
  WRITE (6, '(3x,a,i4,a)')    'hydradj   = ', hydradj,  ','
  WRITE (6, '(3x,a,i4,a)')    'wndadj    = ', wndadj,   ','
  WRITE (6, '(3x,a,i4,a)')    'obropt    = ', obropt,   ','
  WRITE (6, '(3x,a,f16.4,a)') 'obrzero   = ', obrzero,  ','
  WRITE (6, '(3x,a,i4,a)')    'ext_lbc   = ', ext_lbc,  ','
  WRITE (6, '(3x,a,i4,a)')    'ext_vbc   = ', ext_vbc,  ','
  WRITE (6, '(3x,a,i4,a)')    'exttrnopt = ', exttrnopt,','
!wdt-williams 2002-01-17 GMB ntmerge
  WRITE (6, '(3x,a,i4,a)')    'extntmrg  = ', extntmrg, ','
  WRITE (6, '(1x,a)')      '&END'

!  2  CONTINUE

!-----------------------------------------------------------------------
!
!  Allocate and initalize arrays based on dimension parameters
!  read in from the input file
!
!-----------------------------------------------------------------------

  select case ( extdopt )

  case ( 0 )  ! ARPS grid

  select case ( extdfmt )
  case ( 1 )
  FMT = 'bin'
  case ( 2 )
  FMT = 'asc'
  case ( 3 )
  FMT = 'hdf'
  case ( 4 )
  FMT = 'pak'
  case ( 6 )
  FMT = 'bn2'
  case ( 7 )
  FMT = 'net'
  case ( 8 )
  FMT = 'npk'
  case ( 10 )
  FMT = 'grb'
  case default
  PRINT*,'Invalid input history data format given by extdfmt.'
  STOP
  END select

  tmp_ch=trim(dir_extd)//trim(extdname)//'.'//FMT//'grdbas'
  PRINT*,'tmp_ch =', trim(tmp_ch)

  CALL get_dims_from_data(extdfmt,                                      &
       trim(dir_extd)//trim(extdname)//'.'//FMT//'grdbas',              &
       nx_ext,ny_ext,nz_ext,nstyp_ext, ireturn)
  nstyp_ext = min(nstyp_ext,nstyps)

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

  case ( 1 )  ! RUC (Hybrid-B, GRIB #87) from NMC is 81x62x25
  nx_ext=81
  ny_ext=62
  nz_ext=25

  case ( 2 )  ! ETA (GRIB #212, 40km) from NMC is 185x129x39
  nx_ext = 185
  ny_ext = 129
  nz_ext = 39

  case ( 3 ) ! OLAPS for 95 is 91x73x41
  nx_ext=91
  ny_ext=73
  nz_ext=41

  case ( 4 ) ! RUC in GEMPAK format from Rossby (e.g., those of 1995)
  nx_ext=81
  ny_ext=62
  nz_ext=41

  case ( 5 ) ! GEMPAK ETA data

  PRINT*,'Need to find out the grid size and set in ext2arps.f.'
  STOP

!    nx_ext = ??
!    ny_ext = ??
!    nz_ext = ??

  case ( 6 ) ! Special COAMPS files

  PRINT*,'Please specify the size of COAMPS input data grid '
  PRINT*,'by editing file ext2arps.f.'
  STOP

!    nx_ext = ??
!    ny_ext = ??
!    nz_ext = ??

  case ( 7 ) ! Isobaric RUC2 GRIB file (AWIPS #211) from OSO
  nx_ext=93
  ny_ext=65
  nz_ext=37

  case ( 8 ) ! Global Reanalysis on T62 Gaussian grid
  nx_ext = 192
  ny_ext = 94
  nz_ext = 28

  case ( 9 ) ! RUC-2 (GEMPAK) from Doplight is 151x113x41
  nx_ext=151
  ny_ext=113
  nz_ext=41

  case ( 10 )  ! ETA (GEMPAK #104, 80km) from NMC is 147x110x39
  nx_ext=147
  ny_ext=110
  nz_ext=39

  case ( 11 ) ! Native RUC2 GRIB file (Grid #236) from OSO
  nx_ext=151
  ny_ext=113
  nz_ext=40

  case ( 12 ) ! Isobaric RUC2 GRIB file (Grid #236) from OSO
  nx_ext=151
  ny_ext=113
  nz_ext=37

  case ( 13 )  ! AVN (GRIB #3, 1x1 degree) from NCEP is 360x181x26
  nx_ext = 360
  ny_ext = 181
  nz_ext = 26

  case ( 14 ) ! Binary format lat/lon data
  nx_ext = 26
  ny_ext = 16
  nz_ext = 26

  case default

  WRITE(6,'(/a/)')                                                      &
      'extdopt option was invalid. Please specify a new option.'
  STOP

  END select

  PRINT*,'nx_ext, ny_ext, nz_ext = ', nx_ext, ny_ext, nz_ext

  allocate(  allocate(  allocate(
  allocate(xscl(nx),stat=istatus)
  allocate(yscl(ny),stat=istatus)
  allocate(zp(nx,ny,nz),stat=istatus)

  allocate(zs(nx,ny,nz),stat=istatus)
  allocate(j1(nx,ny,nz),stat=istatus)
  allocate(j2(nx,ny,nz),stat=istatus)
  allocate(j3(nx,ny,nz),stat=istatus)
  allocate(aj3z(nx,ny,nz),stat=istatus)

  allocate(hterain(nx,ny),stat=istatus)
!wdt update
  allocate(htrn_ext(nx,ny),stat=istatus)
  allocate(mapfct(nx,ny,8),stat=istatus)

  allocate(  allocate(  allocate(  allocate(pprt(nx,ny,nz),stat=istatus)
  allocate(ptprt(nx,ny,nz),stat=istatus)
  allocate(qv(nx,ny,nz),stat=istatus)
  allocate(qc(nx,ny,nz),stat=istatus)
  allocate(qr(nx,ny,nz),stat=istatus)
  allocate(qi(nx,ny,nz),stat=istatus)
  allocate(qs(nx,ny,nz),stat=istatus)
  allocate(qh(nx,ny,nz),stat=istatus)
!
  allocate(pbar(nx,ny,nz),stat=istatus)
  allocate(ptbar(nx,ny,nz),stat=istatus)
  allocate(qvbar(nx,ny,nz),stat=istatus)
  allocate(ubar(nx,ny,nz),stat=istatus)
  allocate(vbar(nx,ny,nz),stat=istatus)
  allocate(wbar(nx,ny,nz),stat=istatus)
  allocate(rhobar(nx,ny,nz),stat=istatus)
  allocate(rhostr(nx,ny,nz),stat=istatus)
  allocate(wcont(nx,ny,nz),stat=istatus)
  allocate(csndsq(nx,ny,nz),stat=istatus)

  allocate(tsfc(nx,ny,0:nstyps),stat=istatus)
  allocate(tsoil(nx,ny,0:nstyps),stat=istatus)
  allocate(wetsfc(nx,ny,0:nstyps),stat=istatus)
  allocate(wetdp(nx,ny,0:nstyps),stat=istatus)
  allocate(wetcanp(nx,ny,0:nstyps),stat=istatus)
  allocate(snowdpth(nx,ny),stat=istatus)

  allocate(soiltyp (nx,ny,nstyps),stat=istatus)
  allocate(stypfrct(nx,ny,nstyps),stat=istatus)
  allocate(vegtyp (nx,ny),stat=istatus)
  allocate(lai    (nx,ny),stat=istatus)
  allocate(roufns (nx,ny),stat=istatus)
  allocate(veg    (nx,ny),stat=istatus)

  allocate(tem1(nx,ny,nz),stat=istatus)
  allocate(tem2(nx,ny,nz),stat=istatus)
  allocate(tem3(nx,ny,nz),stat=istatus)
  allocate(tem4(nx,ny,nz),stat=istatus)
  allocate(tem5(nx,ny,nz),stat=istatus)
  allocate(tem6(nx,ny,nz),stat=istatus)
  allocate(tem7(nx,ny,nz),stat=istatus)
  allocate(tem8(nx,ny,nz),stat=istatus)

  allocate(xs2d(nx,ny),stat=istatus)
  allocate(ys2d(nx,ny),stat=istatus)
  allocate(xu2d(nx,ny),stat=istatus)
  allocate(yu2d(nx,ny),stat=istatus)
  allocate(xv2d(nx,ny),stat=istatus)
  allocate(yv2d(nx,ny),stat=istatus)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: xw & yw
  allocate(xw(nx,ny),stat=istatus)
  allocate(yw(nx,ny),stat=istatus)

  allocate(iscl(nx,ny),stat=istatus)
  allocate(jscl(nx,ny),stat=istatus)
  allocate(iu(nx,ny),stat=istatus)
  allocate(ju(nx,ny),stat=istatus)
  allocate(iv(nx,ny),stat=istatus)
  allocate(jv(nx,ny),stat=istatus)


  x=0.0
  y=0.0
  z=0.0

  xscl=0.0
  yscl=0.0
  zp=0.0

  zs=0.0
  j1=0.0
  j2=0.0
  j3=0.0
  aj3z=0.0
  hterain=0.0
!wdt update
  htrn_ext=0.0
  mapfct=0.0

  u=0.0
  v=0.0
  w=0.0
  pprt=0.0
  ptprt=0.0
  qv=0.0
  qc=0.0
  qr=0.0
  qi=0.0
  qs=0.0
  qh=0.0
!
  pbar=0.0
  ptbar=0.0
  qvbar=0.0
  ubar=0.0
  vbar=0.0
  wbar=0.0
  rhobar=0.0
  rhostr=0.0
  wcont=0.0
  csndsq=0.0

  tsfc=0.0
  tsoil=0.0
  wetsfc=0.0
  wetdp=0.0
  wetcanp=0.0
  snowdpth=0.0

  soiltyp =0
  stypfrct=0.0
  stypfrct(:,:,1) =1.0
  vegtyp =0
  lai    =0.0
  roufns =0.0
  veg    =0.0

  tem1=0.0
  tem2=0.0
  tem3=0.0
  tem4=0.0
  tem5=0.0
  tem6=0.0
  tem7=0.0
  tem8=0.0

  xs2d=0.0
  ys2d=0.0
  xu2d=0.0
  yu2d=0.0
  xv2d=0.0
  yv2d=0.0

  iscl=0
  jscl=0
  iu=0
  ju=0
  iv=0
  jv=0
!
!-----------------------------------------------------------------------
!
!  Allocate and initialize external grid variables
!
!-----------------------------------------------------------------------
!
  allocate(x_ext(nx_ext),stat=istatus)
  allocate(y_ext(ny_ext),stat=istatus)
  allocate(xu_ext(nx_ext),stat=istatus)
  allocate(yu_ext(ny_ext),stat=istatus)
  allocate(xv_ext(nx_ext),stat=istatus)
  allocate(yv_ext(ny_ext),stat=istatus)
  allocate(lat_ext(nx_ext,ny_ext),stat=istatus)
  allocate(lon_ext(nx_ext,ny_ext),stat=istatus)
  allocate(latu_ext(nx_ext,ny_ext),stat=istatus)
  allocate(lonu_ext(nx_ext,ny_ext),stat=istatus)
  allocate(latv_ext(nx_ext,ny_ext),stat=istatus)
  allocate(lonv_ext(nx_ext,ny_ext),stat=istatus)
  allocate(zs_ext(nx,ny,nz_ext),stat=istatus)
!
  allocate(p_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(hgt_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(zp2_ext(nx,ny,nz_ext),stat=istatus)
  allocate(zp_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(t_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(u_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(v_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(w_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qv_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qc_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qr_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qi_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qs_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(qh_ext(nx_ext,ny_ext,nz_ext),stat=istatus)

  allocate(tsfc_ext   (nx_ext,ny_ext,0:nstyps),stat=istatus)
  allocate(tsoil_ext  (nx_ext,ny_ext,0:nstyps),stat=istatus)
  allocate(wetsfc_ext (nx_ext,ny_ext,0:nstyps),stat=istatus)
  allocate(wetdp_ext  (nx_ext,ny_ext,0:nstyps),stat=istatus)
  allocate(wetcanp_ext(nx_ext,ny_ext,0:nstyps),stat=istatus)
  allocate(snowdpth_ext(nx_ext,ny_ext),stat=istatus)

  allocate(soiltyp_ext (nx_ext,ny_ext,nstyps),stat=istatus)
  soiltyp_ext = 0
  allocate(stypfrct_ext(nx_ext,ny_ext,nstyps),stat=istatus)
  allocate(vegtyp_ext (nx_ext,ny_ext),stat=istatus)
  allocate(lai_ext    (nx_ext,ny_ext),stat=istatus)
  allocate(roufns_ext (nx_ext,ny_ext),stat=istatus)
  allocate(veg_ext    (nx_ext,ny_ext),stat=istatus)

  allocate(trn_ext    (nx_ext,ny_ext),stat=istatus)
  allocate(psfc_ext   (nx_ext,ny_ext),stat=istatus)

  allocate(t_2m_ext (nx_ext,ny_ext),stat=istatus)
  allocate(rh_2m_ext(nx_ext,ny_ext),stat=istatus)
  allocate(u_10m_ext(nx_ext,ny_ext),stat=istatus)
  allocate(v_10m_ext(nx_ext,ny_ext),stat=istatus)

  allocate(dxfld(nx_ext),stat=istatus)
  allocate(dyfld(ny_ext),stat=istatus)
  allocate(rdxfld(nx_ext),stat=istatus)
  allocate(rdyfld(ny_ext),stat=istatus)
  allocate(dxfldu(nx_ext),stat=istatus)
  allocate(dyfldu(ny_ext),stat=istatus)
  allocate(rdxfldu(nx_ext),stat=istatus)
  allocate(rdyfldu(ny_ext),stat=istatus)
  allocate(dxfldv(nx_ext),stat=istatus)
  allocate(dyfldv(ny_ext),stat=istatus)
  allocate(rdxfldv(nx_ext),stat=istatus)
  allocate(rdyfldv(ny_ext),stat=istatus)
  allocate(tem1_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(tem2_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(tem3_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(tem4_ext(nx_ext,ny_ext,nz_ext),stat=istatus)
  allocate(tem5_ext(nx_ext,ny_ext,nz_ext),stat=istatus)

  allocate(xa_ext(nx_ext,ny_ext),stat=istatus)
  allocate(ya_ext(nx_ext,ny_ext),stat=istatus)
  allocate(avgzs_ext(nx,ny,nz_ext),stat=istatus)

!
  x_ext=0.0
  y_ext=0.0
  xu_ext=0.0
  yu_ext=0.0
  xv_ext=0.0
  yv_ext=0.0
  lat_ext=0.0
  lon_ext=0.0
  latu_ext=0.0
  lonu_ext=0.0
  latv_ext=0.0
  lonv_ext=0.0
  zs_ext=0.0
!
  p_ext=0.0
  hgt_ext=0.0
  zp2_ext=0.0
  zp_ext=0.0
  t_ext=0.0
  u_ext=0.0
  v_ext=0.0
  w_ext=0.0
  qv_ext=0.0
  qc_ext=0.0
  qr_ext=0.0
  qi_ext=0.0
  qs_ext=0.0
  qh_ext=0.0

  tsfc_ext    =-999.0
  tsoil_ext   =-999.0
  wetsfc_ext  =-999.0
  wetdp_ext   =-999.0
  wetcanp_ext =-999.0
  snowdpth_ext=-999.0

  soiltyp = -999.0
  stypfrct= -999.0
  vegtyp = -999
  lai    = -999.0
  roufns = -999.0
  veg    = -999.0

  trn_ext    =0.0
  psfc_ext   =0.0

  t_2m_ext =0.0
  rh_2m_ext=0.0
  u_10m_ext=0.0
  v_10m_ext=0.0

  dxfld=0.0
  dyfld=0.0
  rdxfld=0.0
  rdyfld=0.0
  dxfldu=0.0
  dyfldu=0.0
  rdxfldu=0.0
  rdyfldu=0.0
  dxfldv=0.0
  dyfldv=0.0
  rdxfldv=0.0
  rdyfldv=0.0
  tem1_ext=0.0
  tem2_ext=0.0
  tem3_ext=0.0
  tem4_ext=0.0
  tem5_ext=0.0

  xa_ext=0.0
  ya_ext=0.0
  avgzs_ext=0.0
!
!-----------------------------------------------------------------------
!
!  Set up ARPS map projection
!
!-----------------------------------------------------------------------
!
  latnot(1)=trulat1
  latnot(2)=trulat2
  CALL setmapr(mapproj,sclfct,latnot,trulon)
!
!-----------------------------------------------------------------------
!
!  Set up ARPS grid
!
!-----------------------------------------------------------------------
!
  IF (exttrnopt == 1) THEN ! don't really do grid here, just set map proj
    ternopt = -1
    hterain = 0
  END IF

  CALL inigrd(nx,ny,nz,                                                 &
              x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)

  DO k=2,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Set up the height levels on which to define the base-state.
!
!-----------------------------------------------------------------------
!
  deltaz = depthp/(lvlprof-1)

  WRITE(6,'(/a,/a,f7.2,a/)')                                            &
      ' The base state is formed from a mean sounding created',         &
      ' with a ',deltaz,' meter interval.'

  DO k=1,lvlprof
    zsnd(k) = (k-1)*deltaz
  END DO
!
!-----------------------------------------------------------------------
!
!  Loop through the data times provided via NAMELIST.
!
!-----------------------------------------------------------------------
!
  iniotfu = 21  ! FORTRAN unit number used for data output

  DO ifile=1,nextdfil
!
!-----------------------------------------------------------------------
!
!  Time conversions.
!  Formats:  extdtime='1994-05-06.18:00:00+000:00:00'
!            julfname='941261800'
!
!-----------------------------------------------------------------------
!
    READ(extdtime(ifile),'(a19,1x,a9)') extdinit,extdfcst
    IF(extdfcst == '         ') extdfcst='000:00:00'
    READ(extdinit,                                                      &
        '(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)',ERR=920,END=920)           &
        iyr,imo,iday,ihr,imin,isec
    CALL julday(iyr,imo,iday,jldy)
    myr=MOD(iyr,100)
    ifhr=0
    ifmin=0
    ifsec=0
    READ(extdfcst,                                                      &
        '(i3,1x,i2,1x,i2)',ERR=4,END=4) ifhr,ifmin,ifsec
    4    CONTINUE
    mfhr=MOD(ifhr,24)
    jldy = jldy + ifhr/24
    WRITE(julfname,                                                     &
         '(i2.2,i3.3,i2.2,i2.2)') myr,jldy,ihr,mfhr
    CALL ctim2abss(iyr,imo,iday,ihr,imin,isec,iabssec)
    jabssec=(ifhr*3600) + (ifmin*60) + ifsec + iabssec
    kftime=jabssec-initsec
    curtim=FLOAT(kftime)
    WRITE(6,'(a,a9,a,/19x,a,a19,a/a,a/a,i16,a,/,i26,a)')                &
        ' Calling rdextfil, looking for ',                              &
         extdfcst,' hour forecast ',                                    &
        'initialized at ',extdinit,' UTC',                              &
        ' Julian filename: ',julfname,                                  &
        ' Which is ',jabssec,' abs seconds or ',kftime,                 &
        ' seconds from the ARPS initial time.'
!
!-----------------------------------------------------------------------
!
!  Call getextd to reads and converts data to ARPS units
!
!-----------------------------------------------------------------------
!

    select case ( extdopt )

    case ( 0 )  ! ARPS grid

      CALL getarps(nx_ext,ny_ext,nz_ext,                                  &
                   dir_extd,extdname,extdopt,extdfmt,                     &
                   extdinit,extdfcst,julfname,nstyps,                     &
                   iproj_ext,scale_ext,                                   &
                   trlon_ext,latnot_ext,x0_ext,y0_ext,                    &
                   lat_ext,lon_ext,latu_ext,lonu_ext,latv_ext,lonv_ext,   &
                   p_ext,hgt_ext,zp_ext,t_ext,qv_ext,  &
                   u_ext,tem4_ext,v_ext,tem3_ext,w_ext,   &
                   qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                    &
                   soiltyp_ext,stypfrct_ext,vegtyp_ext,                   &
                   lai_ext,roufns_ext,veg_ext,                            &
                   tsfc_ext,tsoil_ext,                                    &
                   wetsfc_ext,wetdp_ext,wetcanp_ext,                      &
                   tem1_ext,tem2_ext,istatus)
      ! tem3_ext holds uatv_ext and tem4_ext holds vatu_ext until
      ! u & v rotated to new map projection below (see calls to uvetomp)
      trn_ext(:,:) = zp_ext(:,:,2)

    case ( 1 )
!
!-----------------------------------------------------------------------
!
!  Get NMC RUC GRIB #87
!
!-----------------------------------------------------------------------
!
      CALL getnmcruc87(nx_ext,ny_ext,nz_ext,                              &
                       dir_extd,extdname,extdopt,extdfmt,                 &
                       extdinit,extdfcst,julfname,                        &
                       iproj_ext,scale_ext,                               &
                       trlon_ext,latnot_ext,x0_ext,y0_ext,                &
                       lat_ext,lon_ext,                                   &
                       p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,            &
                       qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                &
                       tsfc_ext,tsoil_ext,                                &
                       wetsfc_ext,wetdp_ext,wetcanp_ext,                  &
                       trn_ext,psfc_ext,                                  &
                       istatus)

    case ( 2 )
!
!-----------------------------------------------------------------------
!
!  Get NMC ETA GRIB #212
!
!-----------------------------------------------------------------------
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. soiltyp_ext
      CALL getnmceta212(nx_ext,ny_ext,nz_ext,                             &
                        dir_extd,extdname,extdopt,extdfmt,                &
                        extdinit,extdfcst,julfname,                       &
                        iproj_ext,scale_ext,                              &
                        trlon_ext,latnot_ext,x0_ext,y0_ext,               &
                        lat_ext,lon_ext,                                  &
                        p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,           &
                        qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,               &
                        tsfc_ext,tsoil_ext,                               &
                        wetsfc_ext,wetdp_ext,wetcanp_ext,                 &
                        snowdpth_ext,trn_ext,psfc_ext,soiltyp_ext,        &
                        istatus)
      nstyp_ext = 1
      stypfrct_ext(:,:,1) = 1.
      stypfrct_ext(:,:,2:nstyps) = 0.


    case ( 3 )
!
!-----------------------------------------------------------------------
!
!  Get LAPS data. Need LAPS library (see makearps for ext2arps.laps)
!
!-----------------------------------------------------------------------
!
      CALL getlaps(nx_ext,ny_ext,nz_ext,dir_extd,                         &
                   extdinit,extdfcst,julfname,iabssec,                    &
                   iproj_ext,scale_ext,                                   &
                   trlon_ext,latnot_ext,x0_ext,y0_ext,                    &
                   lat_ext,lon_ext,                                       &
                   p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                &
                   qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                    &
                   istatus)

    case ( 4 )
!
!-----------------------------------------------------------------------
!
!  Get RUC/MAPS data in GEMPAK format with a special setup for GEMPAK
!  path and library (see makearps for ext2arps.gemruc).
!
!-----------------------------------------------------------------------
!
      CALL getgemruc(nx_ext,ny_ext,nz_ext,dir_extd,extdname,              &
                     extdinit,extdfcst,julfname,iabssec,                  &
                     iproj_ext,scale_ext,                                 &
                     trlon_ext,latnot_ext,x0_ext,y0_ext,                  &
                     lat_ext,lon_ext,                                     &
                     p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,              &
                     qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                  &
                     istatus, tem1_ext)

    case ( 5 )

      CALL getgemeta(nx_ext,ny_ext,nz_ext,dir_extd,extdname,              &
                     extdinit,extdfcst,julfname,iabssec,                  &
                     iproj_ext,scale_ext,                                 &
                     trlon_ext,latnot_ext,x0_ext,y0_ext,                  &
                     lat_ext,lon_ext,                                     &
                     p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,              &
                     qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                  &
                     istatus, tem1_ext)

    case ( 6 )
!
!-----------------------------------------------------------------------
!
!  Get COAMPS data
!
!-----------------------------------------------------------------------
!
      CALL getcoamps(nx_ext, ny_ext, nz_ext, dir_extd,                    &
                     extdinit,extdfcst,                                   &
                     iproj_ext,scale_ext,                                 &
                     trlon_ext,latnot_ext,x0_ext,y0_ext,                  &
                     lat_ext,lon_ext,                                     &
                     p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,              &
                     qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                  &
                     tsfc_ext,tsoil_ext,                                  &
                     wetsfc_ext,wetdp_ext,wetcanp_ext,istatus)

    case ( 7 )
!
!-----------------------------------------------------------------------
!
!  Get NMC RUC2 GRIB #211
!
!-----------------------------------------------------------------------
!
      CALL getnmcruc211(nx_ext,ny_ext,nz_ext,                             &
                        dir_extd,extdname,extdopt,extdfmt,                &
                        extdinit,extdfcst,julfname,                       &
                        iproj_ext,scale_ext,                              &
                        trlon_ext,latnot_ext,x0_ext,y0_ext,               &
                        lat_ext,lon_ext,                                  &
                        p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,           &
                        qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,               &
                        tsfc_ext,tsoil_ext,                               &
                        wetsfc_ext,wetdp_ext,wetcanp_ext,                 &
                        trn_ext,psfc_ext, t_2m_ext, rh_2m_ext,            &
                        u_10m_ext, v_10m_ext, istatus)

    case ( 8 )
!
!-----------------------------------------------------------------------
!
!  Get NCEP global re-analysis on T62 Gaussian lat/lon grid
!
!-----------------------------------------------------------------------
!
      CALL getreanalt62(nx_ext,ny_ext,nz_ext,                             &
                        dir_extd,extdname,extdopt,extdfmt,                &
                        extdinit,extdfcst,julfname,                       &
                        iproj_ext,scale_ext,                              &
                        trlon_ext,latnot_ext,x0_ext,y0_ext,               &
                        lat_ext,lon_ext,                                  &
                        p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,           &
                        qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,               &
                        tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,          &
                        wetcanp_ext, trn_ext,psfc_ext, istatus)

    case ( 9 )
!
!-----------------------------------------------------------------------
!
!  Get NMC RUC2 GEM
!
!-----------------------------------------------------------------------
!
      CALL getgemruc2(nx_ext,ny_ext,nz_ext,dir_extd,extdname,             &
                     extdinit,extdfcst,julfname,iabssec,                  &
                     iproj_ext,scale_ext,                                 &
                     trlon_ext,latnot_ext,x0_ext,y0_ext,                  &
                     lat_ext,lon_ext,                                     &
                     p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,              &
                     qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                  &
                     istatus, tem1_ext)

    case ( 10 )
!
!-----------------------------------------------------------------------
!
!  Get NMC ETA GEMPAK #104
!
!-----------------------------------------------------------------------
!
      CALL getgemeta2(nx_ext,ny_ext,nz_ext,dir_extd,extdname,             &
                      extdinit,extdfcst,julfname,iabssec,                 &
                      iproj_ext,scale_ext,                                &
                      trlon_ext,latnot_ext,x0_ext,y0_ext,                 &
                      lat_ext,lon_ext,                                    &
                      p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,             &
                      qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                 &
                      istatus, tem1_ext)


!
!-----------------------------------------------------------------------
!
!  Get NCEP RUC2 Native Coordinate GRIB (#236)
!
!-----------------------------------------------------------------------
!

    case ( 11 )

      CALL getnmcrucn236(nx_ext,ny_ext,nz_ext,                            &
                         dir_extd,extdname,extdopt,extdfmt,               &
                         extdinit,extdfcst,julfname,                      &
                         iproj_ext,scale_ext,                             &
                         trlon_ext,latnot_ext,x0_ext,y0_ext,              &
                         lat_ext,lon_ext,                                 &
                         p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,          &
                         qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,              &
                         tsfc_ext,tsoil_ext,                              &
                         wetsfc_ext,wetdp_ext,wetcanp_ext,                &
                         trn_ext,psfc_ext,snowdpth_ext,                   &
                         istatus)

!
!-----------------------------------------------------------------------
!
!  Get NCEP RUC2 Isobaric Coordinate GRIB (#236)
!
!-----------------------------------------------------------------------
!

    case ( 12 )

      CALL getnmcrucp236(nx_ext,ny_ext,nz_ext,                            &
                        dir_extd,extdname,extdopt,extdfmt,                &
                        extdinit,extdfcst,julfname,                       &
                        iproj_ext,scale_ext,                              &
                        trlon_ext,latnot_ext,x0_ext,y0_ext,               &
                        lat_ext,lon_ext,                                  &
                        p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,           &
                        qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,               &
                        tsfc_ext,tsoil_ext,                               &
                        wetsfc_ext,wetdp_ext,wetcanp_ext,                 &
                        trn_ext,psfc_ext, t_2m_ext, rh_2m_ext,            &
                        u_10m_ext, v_10m_ext, snowdpth_ext,               &
                        istatus)

    case ( 13 )
!
!-----------------------------------------------------------------------
!
!  Get NCEP AVN GRIB #3
!
!-----------------------------------------------------------------------
!
      CALL getncepavn3(nx_ext,ny_ext,nz_ext,                              &
                       dir_extd,extdname,extdopt,extdfmt,                 &
                       extdinit,extdfcst,julfname,                        &
                       iproj_ext,scale_ext,                               &
                       trlon_ext,latnot_ext,x0_ext,y0_ext,                &
                       lat_ext,lon_ext,                                   &
                       p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,            &
                       qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                &
                       tsfc_ext,tsoil_ext,                                &
                       wetsfc_ext,wetdp_ext,wetcanp_ext,                  &
                       snowdpth_ext,trn_ext,psfc_ext,                     &
                       istatus)

    case ( 14 ) ! Binary format lat/lon data

      CALL get_avn_bin(nx_ext,ny_ext,nz_ext,                              &
           dir_extd,extdname,extdinit,extdfcst,julfname,                  &
           iproj_ext,scale_ext,trlon_ext,latnot_ext,x0_ext,y0_ext,        &
           lat_ext,lon_ext,p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,        &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                            &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,           &
           snowdpth_ext,trn_ext,psfc_ext,                                 &
           istatus)

    case default

      WRITE(6,'(/a/)')  &
        'extdopt option was invalid. Please specify a new option.'
      STOP

    END select

    IF(istatus /= 1) GO TO 999

    PRINT*,' '
    CALL a3dmax0(lat_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'lat_ext_min= ', amin,', lat_ext_max=',amax
    CALL a3dmax0(lon_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'lon_ext_min= ', amin,', lon_ext_max=',amax
!
    CALL a3dmax0(p_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'p_ext_min  = ', amin,', p_ext_max  =',amax
    CALL a3dmax0(hgt_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'hgt_ext_min= ', amin,', hgt_ext_max=',amax
    CALL a3dmax0(t_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          't_ext_min  = ', amin,', t_ext_max  =',amax
    CALL a3dmax0(u_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'u_ext_min  = ', amin,', u_ext_max  =',amax
    CALL a3dmax0(v_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'v_ext_min  = ', amin,', v_ext_max  =',amax
    CALL a3dmax0(w_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'w_ext_min  = ', amin,', w_ext_max  =',amax
    CALL a3dmax0(qv_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qv_ext_min = ', amin,', qv_ext_max =',amax
    CALL a3dmax0(qc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qc_ext_min = ', amin,', qc_ext_max =',amax
    CALL a3dmax0(qr_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qr_ext_min = ', amin,', qr_ext_max =',amax
    CALL a3dmax0(qi_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qi_ext_min = ', amin,', qi_ext_max =',amax
    CALL a3dmax0(qs_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qs_ext_min = ', amin,', qs_ext_max =',amax
    CALL a3dmax0(qh_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,             &
               1,nz_ext,1,nz_ext,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'qh_ext_min = ', amin,', qh_ext_max =',amax

!wdt note: FIXME add soiltype, etc,  plus add nstyp
    CALL a3dmax0(tsfc_ext   ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'tsfc_ext_min= ', amin,', tsfc_ext_max=',amax
    CALL a3dmax0(tsoil_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'tsoil_ext_min= ', amin,', tsoil_ext_max=',amax
    CALL a3dmax0(wetsfc_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'wetsfc_ext_min= ', amin,', wetsfc_ext_max=',amax
    CALL a3dmax0(wetdp_ext  ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'wetdp_ext_min= ', amin,', wetdp_ext_max=',amax
    CALL a3dmax0(wetcanp_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'wetcanp_ext_min= ', amin,', wetcanp_ext_max=',amax
    CALL a3dmax0(snowdpth_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,        &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'snowd_ext_min= ', amin,', snow_ext_max=',amax

    CALL a3dmax0(trn_ext    ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'trn_ext_min= ', amin,', trn_ext_max=',amax
    CALL a3dmax0(psfc_ext   ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,         &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'psfc_ext_min= ', amin,', psfc_ext_max=',amax

    CALL a3dmax0(t_2m_ext ,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,           &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'T_2m_ext_min= ', amin,', T_2m_ext_max=',amax
    CALL a3dmax0(rh_2m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,           &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'RH_2m_ext_min= ', amin,', RH_2m_ext_max=',amax
    CALL a3dmax0(u_10m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,           &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'U_10m_ext_min= ', amin,', U_10m_ext_max=',amax
    CALL a3dmax0(v_10m_ext,1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,           &
               1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
          'V_10m_ext_min= ', amin,', V_10m_ext_max=',amax
    PRINT*,' '
!
!-----------------------------------------------------------------------
!
!    First time through the time loop, calculate grid
!    transformation info.
!
!-----------------------------------------------------------------------
!
    IF(ifile == 1) THEN

      DO i=1,nx-1
        xscl(i)=0.5*(x(i)+x(i+1))
      END DO
      xscl(nx)=2.*xscl(nx-1) - xscl(nx-2)
      DO j=1,ny-1
        yscl(j)=0.5*(y(j)+y(j+1))
      END DO
      yscl(ny)=2.*yscl(ny-1) - yscl(ny-2)
!
!-----------------------------------------------------------------------
!
!  Find lat,lon locations of ARPS scalar, u and v grids.
!
!-----------------------------------------------------------------------
!
      CALL xytoll(nx,ny,xscl,yscl,tem1,tem2)
      CALL xytoll(nx,ny,   x,yscl,tem3,tem4)
      CALL xytoll(nx,ny,xscl,   y,tem5,tem6)
      CALL getmapr(iproj_arps,scale_arps,latnot_arps,                     &
                   trlon_arps,xorig_arps,yorig_arps)
!
!-----------------------------------------------------------------------
!
!  Find x,y locations of external grid.
!
!-----------------------------------------------------------------------
!
      CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
      CALL setorig(1,x0_ext,y0_ext)
      DO j=1,ny_ext
        CALL lltoxy(1,1,lat_ext(1,j),lon_ext(1,j),                        &
                        x_ext(1),y_ext(j))
      END DO
      DO i=1,nx_ext
        CALL lltoxy(1,1,lat_ext(i,1),lon_ext(i,1),                        &
                        x_ext(i),y_ext(1))
      END DO
      IF ( extdopt == 0 ) THEN
        DO j=1,ny_ext
          CALL lltoxy(1,1,latu_ext(1,j),lonu_ext(1,j),                        &
                          xu_ext(1),yu_ext(j))
        END DO
        DO i=1,nx_ext
          CALL lltoxy(1,1,latu_ext(i,1),lonu_ext(i,1),                        &
                          xu_ext(i),yu_ext(1))
        END DO
        DO j=1,ny_ext
          CALL lltoxy(1,1,latv_ext(1,j),lonv_ext(1,j),                        &
                          xv_ext(1),yv_ext(j))
        END DO
        DO i=1,nx_ext
          CALL lltoxy(1,1,latv_ext(i,1),lonv_ext(i,1),                        &
                          xv_ext(i),yv_ext(1))
        END DO
      ENDIF

      IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
        DO i=1,nx_ext
          IF(x_ext(i) < 0.0) x_ext(i)=x_ext(i)+360.
        END DO
      END IF
!
!-----------------------------------------------------------------------
!
!  Find x,y locations of ARPS scalar grid in terms of the
!  external grid.
!
!-----------------------------------------------------------------------
!
      CALL lltoxy(nx,ny,tem1,tem2,xs2d,ys2d)
      IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
        DO i=1,nx
          DO j=1,ny
            IF( xs2d(i,j) < 0. ) xs2d(i,j)= xs2d(i,j) + 360.
          END DO
        END DO
      END IF

      CALL setijloc(nx,ny,nx_ext,ny_ext,xs2d,ys2d,                        &
                    x_ext,y_ext,iscl,jscl)
!
!-----------------------------------------------------------------------
!
!  Find x,y locations of ARPS u grid in terms of the
!  external grid.
!
!-----------------------------------------------------------------------
!
      CALL lltoxy(nx,ny,tem3,tem4,xu2d,yu2d)
      IF( extdopt == 8 .OR. extdopt == 13 ) THEN
        DO i=1,nx
          DO j=1,ny
            IF ( xu2d(i,j) < 0. ) xu2d(i,j)= xu2d(i,j) + 360.
          END DO
        END DO
      END IF
      IF ( extdopt == 0 ) THEN
        CALL setijloc(nx,ny,nx_ext,ny_ext,xu2d,yu2d,                        &
                      xu_ext,yu_ext,iu,ju)
      ELSE
        CALL setijloc(nx,ny,nx_ext,ny_ext,xu2d,yu2d,                        &
                      x_ext,y_ext,iu,ju)
      ENDIF
!
!-----------------------------------------------------------------------
!
!  Find x,y locations of ARPS v grid in terms of the
!  external grid.
!
!-----------------------------------------------------------------------
!
      CALL lltoxy(nx,ny,tem5,tem6,xv2d,yv2d)
      IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
        DO i=1,nx
          DO j=1,ny
            IF ( xv2d(i,j) < 0. ) xv2d(i,j)= xv2d(i,j) + 360.
          END DO
        END DO
      END IF
      IF ( extdopt == 0 ) THEN
        CALL setijloc(nx,ny,nx_ext,ny_ext,xv2d,yv2d,                        &
                      xv_ext,yv_ext,iv,jv)
      ELSE
        CALL setijloc(nx,ny,nx_ext,ny_ext,xv2d,yv2d,                        &
                      x_ext,y_ext,iv,jv)
      ENDIF

      xumin=xu2d(1,1)
      xumax=xu2d(1,1)
      DO j=1,ny-1
        xumin=MIN(xumin,xu2d(1 ,j))
        xumax=MAX(xumax,xu2d(nx,j))
      END DO

      yvmin=yv2d(1,1)
      yvmax=yv2d(1,1)
      DO i=1,nx-1
        yvmin=MIN(yvmin,yv2d(i,1 ))
        yvmax=MAX(yvmax,yv2d(i,ny))
      END DO

      IF(xumin < x_ext(1).OR.xumax > x_ext(nx_ext).OR.                    &
            yvmin < y_ext(1).OR.yvmax > y_ext(ny_ext)) THEN
        WRITE(6,'(a/a)')                                                  &
            ' ARPS domain extends outside available external data',       &
            ' domain, ext2arps aborted.'
        WRITE (6,*) "xu min & ext:",xumin,x_ext(1)
        WRITE (6,*) "xu max & ext:",xumax,x_ext(nx_ext)
        WRITE (6,*) "yv min & ext:",yvmin,y_ext(1)
        WRITE (6,*) "yv max & ext:",yvmax,y_ext(ny_ext)
        STOP
      END IF
!
!-----------------------------------------------------------------------
!
!  Test code, for diagnostic testing.
!  Find x,y of Norman sounding in external grid.
!
!-----------------------------------------------------------------------
!
      CALL lltoxy(1,1,latdiag,londiag,xdiag,ydiag)

      IF ( extdopt == 8 .OR. extdopt == 13 ) THEN
        IF ( xdiag < 0. ) xdiag= xdiag + 360.
      END IF

      dmin=((xdiag-x_ext(1))*(xdiag-x_ext(1))+                            &
            (ydiag-y_ext(1))*(ydiag-y_ext(1)))

      idiag=1
      jdiag=1
!
      DO j=1,ny_ext
        DO i=1,nx_ext
          dd=((xdiag-x_ext(i))*(xdiag-x_ext(i))+                          &
              (ydiag-y_ext(j))*(ydiag-y_ext(j)))
          IF(dd < dmin) THEN
            dmin=dd
            idiag=i
            jdiag=j
          END IF
        END DO
      END DO
      WRITE(6,'(a,f10.4,f10.4,/a,i5,i5,a,f10.4,f10.4)')                   &
            ' Nearest ext pt to diagnostic lat,lon: ',                    &
            latdiag,londiag,                                              &
            ' Diagnostic i,j: ',                                          &
            idiag,jdiag,' lat,lon= ',                                     &
            lat_ext(idiag,jdiag),lon_ext(idiag,jdiag)
      WRITE(6,'(///a,/2x,a)') ' External sounding at idiag,jdiag',        &
          'k   pres   hgt   temp   theta   dewp     u     v     dir    spd'
!
!-----------------------------------------------------------------------
!
!  Convert units of external data and write as a sounding.
!
!-----------------------------------------------------------------------
!
      DO k=nz_ext,1,-1
        pmb=.01*p_ext(idiag,jdiag,k)
        tc=t_ext(idiag,jdiag,k)-273.15
        theta=t_ext(idiag,jdiag,k)*(p0/p_ext(idiag,jdiag,k))**rddcp
        IF( qv_ext(idiag,jdiag,k) > 0.) THEN
          smix=qv_ext(idiag,jdiag,k)/(1.-qv_ext(idiag,jdiag,k))
          e=(pmb*smix)/(0.62197 + smix)
          bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034)
          alge = ALOG(bige/6.112)
          tdc = (alge * 243.5) / (17.67 - alge)
        ELSE
          tdc = tc-30.
        END IF

        CALL uvrotdd(1,1,lon_ext(idiag,jdiag),                            &
                     u_ext(idiag,jdiag,k),v_ext(idiag,jdiag,k),           &
                     dir,spd)

        WRITE(6,'(i4,f6.0,f9.0,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1)')      &
                      k,pmb,                                              &
                      hgt_ext(idiag,jdiag,k),                             &
                      tc,theta,tdc,                                       &
                      u_ext(idiag,jdiag,k),                               &
                      v_ext(idiag,jdiag,k),                               &
                      dir,spd
      END DO
!
!-----------------------------------------------------------------------
!
!  Restore map projection to ARPS grid.
!
!-----------------------------------------------------------------------
!
      CALL setmapr(iproj_arps,scale_arps,latnot_arps,                     &
                   trlon_arps)
      CALL setorig(1,xorig_arps,yorig_arps)
!
!-----------------------------------------------------------------------
!
!  Find the min and max of the external indices that will be used
!  to be sure interpolation is within external dataset.
!
!-----------------------------------------------------------------------
!
      iextmn=iscl(1,1)
      jextmn=jscl(1,1)
      iextmx=iscl(1,1)
      jextmx=jscl(1,1)
      DO j=1,ny
        DO i=1,nx
          iextmn=MIN(iextmn,iscl(i,j))
          jextmn=MIN(jextmn,jscl(i,j))

          iextmx=MAX(iextmx,iscl(i,j))
          jextmx=MAX(jextmx,jscl(i,j))

          iextmn=MIN(iextmn,iu(i,j))
          jextmn=MIN(jextmn,ju(i,j))

          iextmx=MAX(iextmx,iu(i,j))
          jextmx=MAX(jextmx,ju(i,j))

          iextmn=MIN(iextmn,iv(i,j))
          jextmn=MIN(jextmn,jv(i,j))

          iextmx=MAX(iextmx,iv(i,j))
          jextmx=MAX(jextmx,jv(i,j))
        END DO
      END DO

    END IF
!
!-----------------------------------------------------------------------
!
!  External data comes in oriented to true north.
!  Change that to u,v in the new coordinate system.
!  (tem3_ext & tem4_ext from above hold uatv_ext & vatu_ext for use here)
!
!-----------------------------------------------------------------------
!
    IF ( extdopt == 0 ) THEN
      DO k=1,nz_ext
        CALL uvetomp(nx_ext,ny_ext,u_ext(1,1,k),tem4_ext(1,1,k),  &
                     lonu_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
      END DO
      u_ext = tem1_ext
      DO k=1,nz_ext
        CALL uvetomp(nx_ext,ny_ext,tem3_ext(1,1,k),v_ext(1,1,k),  &
                     lonv_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
      END DO
      v_ext = tem2_ext
    ELSE
      DO k=1,nz_ext
        CALL uvetomp(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                     lon_ext,tem1_ext(1,1,k),tem2_ext(1,1,k))
      END DO
      u_ext = tem1_ext
      v_ext = tem2_ext
    ENDIF
!
!-----------------------------------------------------------------------
!
!  Process height data
!
!  Interpolate the external heights horizontally to the
!  ARPS x,y.  This is for scalar pts.
!
!-----------------------------------------------------------------------
!
    CALL setdxdy(nx_ext,ny_ext,                                           &
             1,nx_ext,1,ny_ext,                                           &
             x_ext,y_ext,dxfld,dyfld,rdxfld,rdyfld)
    IF ( extdopt == 0 ) THEN
      CALL setdxdy(nx_ext,ny_ext,                                           &
               1,nx_ext,1,ny_ext,                                           &
               xu_ext,yu_ext,dxfldu,dyfldu,rdxfldu,rdyfldu)
      CALL setdxdy(nx_ext,ny_ext,                                           &
               1,nx_ext,1,ny_ext,                                           &
               xv_ext,yv_ext,dxfldv,dyfldv,rdxfldv,rdyfldv)
      DO k=1,nz_ext
        CALL fldint2d(nx,ny,nx_ext,ny_ext,                                  &
                    1,nx,1,ny,                                              &
                    1,nx_ext,1,ny_ext,                                      &
                    iorder,xs2d,ys2d,zp_ext(1,1,k),                         &
                    x_ext,y_ext,iscl,jscl,                                  &
                    dxfld,dyfld,rdxfld,rdyfld,                              &
                    tem1_ext,tem2_ext,tem3_ext,                             &
                    zp2_ext(1,1,k))
      END DO
    ENDIF
    DO k=1,nz_ext
      CALL fldint2d(nx,ny,nx_ext,ny_ext,                                  &
                  1,nx,1,ny,                                              &
                  1,nx_ext,1,ny_ext,                                      &
                  iorder,xs2d,ys2d,hgt_ext(1,1,k),                        &
                  x_ext,y_ext,iscl,jscl,                                  &
                  dxfld,dyfld,rdxfld,rdyfld,                              &
                  tem1_ext,tem2_ext,tem3_ext,                             &
                  zs_ext(1,1,k))
    END DO
!
!-----------------------------------------------------------------------
!
!  Calculate x and y coordinates of external grid in ARPS coordinate
!  system.  Store them in xa_ext and ya_ext.
!
!-----------------------------------------------------------------------
!
    CALL lltoxy(nx_ext,ny_ext,lat_ext,lon_ext,                            &
                xa_ext,ya_ext)

    IF(ifile == 1) THEN
!
!-----------------------------------------------------------------------
!
!  Interpolate external terrain to arps grid (if exttrnopt=1).
!
!-----------------------------------------------------------------------
!
!wdt update exttrnopt
!wdt-williams 2002-01-17 GMB ntmerge
      IF (exttrnopt > 0) THEN ! set arps terrain to match external terrain

        CALL mkarps2d (nx_ext,ny_ext,nx,ny,  &
                   iorder,iscl,jscl,x_ext,y_ext,  &
                   xs2d,ys2d,trn_ext,htrn_ext,  &
                   dxfld,dyfld,rdxfld,rdyfld,        &
                   tem1_ext(1,1,1),tem1_ext(1,1,2),  &
                   tem1_ext(1,1,3),tem1_ext(1,1,4))

        IF (exttrnopt == 1) THEN
          hterain = htrn_ext
!wdt-williams 2002-01-17 GMB ntmerge
        ELSE 
          ntmergeinv = 1.d0/extntmrg
          DO j = 1,ny
            DO i = 1,nx
              idist = max(0,min(extntmrg,i-2,nx-2-i,j-2,ny-2-j)) 
              mfac = idist*ntmergeinv
              hterain(i,j) = (1.d0-mfac)*htrn_ext(i,j)  &
                              + mfac*hterain(i,j)
            END DO
          END DO
        END IF

        ternopt = -1
        CALL inigrd(nx,ny,nz,  &
                    x,y,z,zp,hterain,mapfct,j1,j2,j3,tem1,tem2,tem3)

        DO k=2,nz-1
          DO j=1,ny-1
            DO i=1,nx-1
              aj3z(i,j,k)=0.5*(j3(i,j,k)+j3(i,j,k-1))
            END DO
          END DO
        END DO

      END IF

!wdt update 2002-01-17
!-----------------------------------------------------------------------
!
!  Write out terrain data
!
!-----------------------------------------------------------------------

      IF (terndmp > 0) THEN

        ternfn = runname(1:lfnkey)//".trndata"
        lternfn = lfnkey + 8

        IF( dirname /= ' ' ) THEN

          temchar = ternfn
          ternfn = dirname(1:ldirnam)//'/'//temchar
          lternfn  = lternfn + ldirnam + 1

        END IF

        CALL fnversn(ternfn, lternfn )

        PRINT *, 'Write terrain data to ',ternfn(1:lternfn)

        CALL writtrn(nx,ny,ternfn(1:lternfn), dx,dy,                &
                   mapproj,trulat1,trulat2,trulon,sclfct,               &
                   ctrlat,ctrlon,hterain)

        !FIXME: add grads control file
        !IF (terndmp == 1)                                               &
        !    CALL trncntl(nx,ny,ternfn(1:lternfn), x1_out,y1_out,   &
        !                 tem11,tem21)

      END IF

!
!-----------------------------------------------------------------------
!
!  Set up z grid at scalar vertical levels.
!
!-----------------------------------------------------------------------
!
      onvf = 0
      CALL avgz(zp, onvf, nx,ny,nz, 1,nx, 1,ny, 1,nz-1, zs)
!
      DO j=1,ny-1
        DO i=1,nx-1
          zs(i,j,nz)=2.*zs(i,j,nz-1) - zs(i,j,nz-2)
        END DO
      END DO

      CALL edgfill(zs,  1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      CALL edgfill(zp,  1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)

    END IF
!
!-----------------------------------------------------------------------
!
!  Data transformations
!  Calculate log of pressure for external grid.
!  Change qv to RHstar   RHStar=sqrt(1.-relative humidity)
!
!-----------------------------------------------------------------------
!
    qvmin=999.
    qvmax=-999.
    DO k=1,nz_ext
      DO j=1,ny_ext
        DO i=1,nx_ext
          qvmax=AMAX1(qvmax,qv_ext(i,j,k))
          qvmin=AMIN1(qvmin,qv_ext(i,j,k))
          tem5_ext(i,j,k)=ALOG(p_ext(i,j,k))
          qvsat=f_qvsat( p_ext(i,j,k), t_ext(i,j,k) )
          qv_ext(i,j,k)=SQRT(AMAX1(0.,(rhmax-(qv_ext(i,j,k)/qvsat))))
!        u_ext(i,j,k)=tem1_ext(i,j,k)
!        v_ext(i,j,k)=tem2_ext(i,j,k)
        END DO
      END DO
    END DO
    PRINT *, ' qv_ext min = ',qvmin
    PRINT *, ' qv_ext max = ',qvmax
!
!-----------------------------------------------------------------------
!
!  Calculate base-state sounding (vertical profile)
!
!-----------------------------------------------------------------------
!
  CALL extmnsnd(nx,ny,nx_ext,ny_ext,nz_ext,lvlprof,                     &
                    iextmn,iextmx,jextmn,jextmx,1,nz_ext,               &
                    xscl,yscl,xa_ext,ya_ext,                            &
                    hgt_ext,tem5_ext,t_ext,                             &
                    qv_ext,u_ext,v_ext,                                 &
                    zsnd,plsnd,psnd,tsnd,ptsnd,rhssnd,qvsnd,            &
                    rhosnd,usnd,vsnd)
!
!-----------------------------------------------------------------------
!
!  Process Pressure data
!
!-----------------------------------------------------------------------
!
    iprtopt=1
    CALL mkarpsvlz(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                   iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,          &
                   hgt_ext,zs_ext,xs2d,ys2d,zs,p_ext,                     &
                   zsnd,plsnd,pbar,pprt,                                  &
                   dxfld,dyfld,rdxfld,rdyfld,                             &
                   tem1_ext,tem2_ext,tem3_ext,                            &
                   tem4_ext)

    DO ksmth=1,nsmooth
      DO k=1,nz
        CALL smooth9p(pprt(1,1,k), nx,ny, 1,nx,1,ny, tem1)
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Process potential temperature data
!
!-----------------------------------------------------------------------
!
    DO k=1,nz_ext
      DO j=1,ny_ext
        DO i=1,nx_ext
          tem5_ext(i,j,k)=t_ext(i,j,k)*((p0/p_ext(i,j,k))**rddcp)
        END DO
      END DO
    END DO
    iprtopt=1
    CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                   iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,          &
                   hgt_ext,zs_ext,xs2d,ys2d,zs,tem5_ext,                  &
                   zsnd,ptsnd,ptbar,ptprt,                                &
                   dxfld,dyfld,rdxfld,rdyfld,                             &
                   tem1_ext,tem2_ext,tem3_ext,                            &
                   tem4_ext)

    DO ksmth=1,nsmooth
      CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                    smfct1,zs,ptprt,tem1,ptprt)
!       CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
!  :                  smfct2,zs,ptprt,tem1,ptprt)
    END DO
!
!-----------------------------------------------------------------------
!
!  Process qv data.
!
!-----------------------------------------------------------------------
!
    iprtopt=0
    CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                   iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,          &
                   hgt_ext,zs_ext,xs2d,ys2d,zs,qv_ext,                    &
                   zsnd,rhssnd,qvbar,qv,                                  &
                   dxfld,dyfld,rdxfld,rdyfld,                             &
                   tem1_ext,tem2_ext,tem3_ext,                            &
                   tem4_ext)

    DO ksmth=1,nsmooth
      CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                    smfct1,zs,   qv,tem1,   qv)
!       CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
!  :                  smfct2,zs,   qv,tem1,   qv)
    END DO
!
!-----------------------------------------------------------------------
!
!  Convert rhstar back to qv for writing, further calculations
!
!-----------------------------------------------------------------------
!
    qvmax=-999.
    qvmin=999.
    DO k=1,nz_ext
      DO j=1,ny_ext
        DO i=1,nx_ext
          qvsat=f_qvsat( p_ext(i,j,k), t_ext(i,j,k) )
          rh=AMAX1(0.01,rhmax-(qv_ext(i,j,k)*qv_ext(i,j,k)))
          qv_ext(i,j,k)=rh*qvsat
          qvmax=AMAX1(qvmax,qv_ext(i,j,k))
          qvmin=AMIN1(qvmin,qv_ext(i,j,k))
        END DO
      END DO
    END DO
    PRINT *, ' qv_ext min = ',qvmin
    PRINT *, ' qv_ext max = ',qvmax
!
    qvmax=-999.
    qvmin=999.
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          pres = pbar(i,j,k)+pprt(i,j,k)
          temp = (ptbar(i,j,k)+ptprt(i,j,k))*((pres/p0)**rddcp)
          qvsat=f_qvsat( pres, temp )
          IF(qvsat > 1.) THEN
            PRINT *, ' qvsat trouble at: ',i,j,k
            PRINT *, ' temp,press,qvsat: ',temp,pres,qvsat
          END IF
          rh=AMAX1(0.01,rhmax-(qv(i,j,k)*qv(i,j,k)))
          rh=AMIN1(rh,rhmax)
          qv(i,j,k)=rh*qvsat
          IF(qv(i,j,k) > 0.5) THEN
            PRINT *, ' qvsat trouble at: ',i,j,k
            PRINT *, ' temp,press,qvsat: ',temp,pres,qv(i,j,k)
            PRINT *, ' rh= ',rh
          END IF
          qvmax=AMAX1(qvmax,qv(i,j,k))
          qvmin=AMIN1(qvmin,qv(i,j,k))
          pres = pbar(i,j,k)
          temp = ptbar(i,j,k)*((pres/p0)**rddcp)
          qvsat=f_qvsat( pres, temp )
          rh=AMAX1(0.01,rhmax-(qvbar(i,j,k)*qvbar(i,j,k)))
          rh=AMIN1(rh,rhmax)
          qvbar(i,j,k)=rh*qvsat
        END DO
      END DO
    END DO
    PRINT *, ' qv     min = ',qvmin
    PRINT *, ' qv     max = ',qvmax
!
!-----------------------------------------------------------------------
!
!  Process qc data.
!  If first data point is less than or equal to -1, then
!  no qc info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
    IF(qc_ext(1,1,1) > -1.0) THEN
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,               &
                     iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,        &
                     hgt_ext,zs_ext,xs2d,ys2d,zs,qc_ext,                  &
                     zsnd,dumsnd,tem1,qc,                                 &
                     dxfld,dyfld,rdxfld,rdyfld,                           &
                     tem1_ext,tem2_ext,tem3_ext,                          &
                     tem4_ext)

      DO ksmth=1,nsmooth
        CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                           &
                    smfct1,zs,   qc,tem1,   qc)
      END DO

    ELSE
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qc(i,j,k)=0.
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Process qr data.
!  If first data point is less than or equal to -1, then
!  no qr info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
    IF(qr_ext(1,1,1) > -1.0) THEN
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,               &
                     iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,        &
                     hgt_ext,zs_ext,xs2d,ys2d,zs,qr_ext,                  &
                     zsnd,dumsnd,tem1,qr,                                 &
                     dxfld,dyfld,rdxfld,rdyfld,                           &
                     tem1_ext,tem2_ext,tem3_ext,                          &
                     tem4_ext)

      DO ksmth=1,nsmooth
        CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                           &
                      smfct1,zs,   qr,tem1,   qr)
      END DO

      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qr(i,j,k)=MAX(qr(i,j,k),0.)
          END DO
        END DO
      END DO

    ELSE
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qr(i,j,k)=0.
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Process qi data.
!  If first data point is less than or equal to -1, then
!  no qi info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
    iceout=0
    IF(qi_ext(1,1,1) > -1.0) THEN
      iceout=1
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,               &
                     iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,        &
                     hgt_ext,zs_ext,xs2d,ys2d,zs,qi_ext,                  &
                     zsnd,dumsnd,tem1,qi,                                 &
                     dxfld,dyfld,rdxfld,rdyfld,                           &
                     tem1_ext,tem2_ext,tem3_ext,                          &
                     tem4_ext)

      DO ksmth=1,nsmooth
        CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                           &
                    smfct1,zs,   qi,tem1,   qi)
      END DO

      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qi(i,j,k)=MAX(qi(i,j,k),0.)
          END DO
        END DO
      END DO

    ELSE
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qi(i,j,k)=0.
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Process qs data.
!  If first data point is less than or equal to -1, then
!  no qs info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
    IF(qs_ext(1,1,1) > -1.0) THEN
      iceout=1
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,               &
                     iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,        &
                     hgt_ext,zs_ext,xs2d,ys2d,zs,qs_ext,                  &
                     zsnd,dumsnd,tem1,qs,                                 &
                     dxfld,dyfld,rdxfld,rdyfld,                           &
                     tem1_ext,tem2_ext,tem3_ext,                          &
                     tem4_ext)

      DO ksmth=1,nsmooth
        CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                           &
                      smfct1,zs,   qs,tem1,   qs)
      END DO

      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qs(i,j,k)=MAX(qs(i,j,k),0.)
          END DO
        END DO
      END DO

    ELSE
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qs(i,j,k)=0.
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Process qh data.
!  If first data point is less than or equal to -1, then
!  no qh info is provided...set to zero.
!
!-----------------------------------------------------------------------
!
    IF(qh_ext(1,1,1) > -1.0) THEN
      iceout=1
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,               &
                     iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,        &
                     hgt_ext,zs_ext,xs2d,ys2d,zs,qh_ext,                  &
                     zsnd,dumsnd,tem1,qh,                                 &
                     dxfld,dyfld,rdxfld,rdyfld,                           &
                     tem1_ext,tem2_ext,tem3_ext,                          &
                     tem4_ext)

      DO ksmth=1,nsmooth
        CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                           &
                      smfct1,zs,   qh,tem1,   qh)
      END DO

      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qh(i,j,k)=MAX(qh(i,j,k),0.)
          END DO
        END DO
      END DO

    ELSE
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qh(i,j,k)=0.
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Process density which has been stuffed into tem2_ext array.
!  We really only need rhobar, so set iprtopt=1 and pass
!  a temporary array in place of rhoprt.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz_ext
      DO j=1,ny_ext
        DO i=1,nx_ext
          tv_ext=t_ext(i,j,k) /                                           &
                 ((1.0+qv_ext(i,j,k))*                                    &
                  (1.0-(qv_ext(i,j,k)/(rddrv+qv_ext(i,j,k)))))
          tem5_ext(i,j,k)=p_ext(i,j,k)/(rd*tv_ext)
        END DO
      END DO
    END DO
!
    iprtopt=1
    CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                   iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,          &
                   hgt_ext,zs_ext,xs2d,ys2d,zs,tem5_ext,                  &
                   zsnd,rhosnd,rhobar,tem1,                               &
                   dxfld,dyfld,rdxfld,rdyfld,                             &
                   tem1_ext,tem2_ext,tem3_ext,                            &
                   tem4_ext)

    DO ksmth=1,nsmooth
      CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                    smfct1,zs,rhobar,tem1,rhobar)
!       CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
!  :                  smfct2,zs,rhobar,tem1,rhobar)
    END DO
!
!-----------------------------------------------------------------------
!
!  Calculate and store the sound wave speed squared in csndsq.
!
!-----------------------------------------------------------------------
!
    IF(csopt == 1) THEN       ! Original definition of sound speed.
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            csndsq(i,j,k)= cpdcv*pbar(i,j,k)/rhobar(i,j,k)
          END DO
        END DO
      END DO

    ELSE IF(csopt == 2) THEN   ! Original sound speed multiplied
                               ! by a factor
      csconst = csfactr**2*cpdcv
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            csndsq(i,j,k)= csconst * pbar(i,j,k)/rhobar(i,j,k)
          END DO
        END DO
      END DO
    ELSE                      ! Specified constant sound speed.
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            csndsq(i,j,k)= csound * csound
          END DO
        END DO
      END DO
    END IF
!
    IF(hydradj == 1 .OR. hydradj == 2) THEN
      pconst=0.5*g*dz
!
!-----------------------------------------------------------------------
!
!  Create thermal buoyancy at each scalar point,
!  which is stored in tem2
!
!-----------------------------------------------------------------------
!
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            qvprt=qv(i,j,k)-qvbar(i,j,k)
            qtot=qc(i,j,k)+qr(i,j,k)+qi(i,j,k)+qs(i,j,k)+qh(i,j,k)
            tem2(i,j,k)=j3(i,j,k)*rhobar(i,j,k)*                          &
                      g*( (ptprt(i,j,k)/ptbar(i,j,k)) +                   &
                          (qvprt/(rddrv+qvbar(i,j,k))) -                  &
                          ((qvprt+qtot)/(1.+qvbar(i,j,k))) )
          END DO
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!  Average thermal buoyancy to w points
!  which is stored in tem3
!
!-----------------------------------------------------------------------
!
      CALL avgsw(tem2,nx,ny,nz,1,nx,1,ny, tem3)

      IF(hydradj == 1) THEN

        DO i=1,nx
          DO j=1,ny
            tem1(i,j,1)=pprt(i,j,1)
            DO k=2,nz-2
              tem1(i,j,k)=                                                &
                    ( (1. - (pconst*j3(i,j,k-1)/csndsq(i,j,k-1)) )*       &
                        tem1(i,j,k-1) +                                   &
                        dz*tem3(i,j,k) ) /                                &
                        (1. + (pconst*j3(i,j,k)/csndsq(i,j,k)) )
              IF(j == 26 .AND. MOD(i,10) == 0) THEN
                IF(k == nz-2) PRINT *, '            Point i= ',i, '  j=26'
                PRINT *, ' k,pprt,tem1,diff =',                           &
                    k,pprt(i,j,k),tem1(i,j,k),                            &
                    (tem1(i,j,k)-pprt(i,j,k))
              END IF
              pprt(i,j,k)=tem1(i,j,k)
            END DO
            pprt(i,j,nz-1)=pprt(i,j,nz-2)
          END DO
        END DO

      ELSE IF(hydradj == 2) THEN

        DO i=1,nx
          DO j=1,ny
            tem1(i,j,nz-1)=pprt(i,j,nz-1)
            DO k=nz-2,2,-1
              tem1(i,j,k)=                                                &
                    ( (1.+ (pconst*j3(i,j,k+1)/csndsq(i,j,k+1)) )*        &
                        tem1(i,j,k+1) -                                   &
                        dz*tem3(i,j,k+1) ) /                              &
                        (1.- (pconst*j3(i,j,k  )/csndsq(i,j,k  )) )
              IF(j == 26 .AND. MOD(i,10) == 0) THEN
                IF(k == nz-2) PRINT *, '            Point i= ',i, '  j=26'
                PRINT *, ' k,pprt,tem1,diff =',                           &
                    k,pprt(i,j,k),tem1(i,j,k),                            &
                    (tem1(i,j,k)-pprt(i,j,k))
              END IF
              pprt(i,j,k)=tem1(i,j,k)
            END DO
            pprt(i,j,1)=pprt(i,j,2)
          END DO
        END DO
      END IF
    ELSE IF (hydradj == 3) THEN
!
!-----------------------------------------------------------------------
!
!  Calculate total pressure, store in tem1.
!  Calculate virtual temperature, store in tem2.
!
!-----------------------------------------------------------------------
!
      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            tem1(i,j,k) = pbar(i,j,k)+pprt(i,j,k)
            temp = (ptbar(i,j,k)+ptprt(i,j,k))*                           &
                   ((tem1(i,j,k)/p0)**rddcp)
            tem2(i,j,k) = temp*(1.0+rvdrd*qv(i,j,k))/(1.0+qv(i,j,k))
          END DO
        END DO
      END DO
!
!-----------------------------------------------------------------------
!
!  Integrate hydrostatic equation, begining at k=2
!
!-----------------------------------------------------------------------
!
      pconst=-g/rd
      DO k=3,nz-1
        DO j=1,ny
          DO i=1,nx
            tvbar=0.5*(tem2(i,j,k)+tem2(i,j,k-1))
            tem1(i,j,k)=tem1(i,j,k-1)*                                    &
                       EXP(pconst*(zs(i,j,k)-zs(i,j,k-1))/tvbar)
            pprt(i,j,k)=tem1(i,j,k)-pbar(i,j,k)
          END DO
        END DO
      END DO
      DO j=1,ny
        DO i=1,nx
          tvbar=0.5*(tem2(i,j,2)+tem2(i,j,1))
          tem1(i,j,1)=tem1(i,j,2)*                                        &
                     EXP(pconst*(zs(i,j,1)-zs(i,j,2))/tvbar)
          pprt(i,j,1)=tem1(i,j,1)-pbar(i,j,1)
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!    Process wind data
!
!-----------------------------------------------------------------------
!
    CALL avgsu(zs_ext,nx,ny,nz_ext, 1,ny, 1,nz_ext, avgzs_ext, tem3_ext)
    CALL avgsu(zs,nx,ny,nz, 1,ny, 1,nz, tem2, tem3)

    DO k=1,nz_ext
      DO j=1,ny
        avgzs_ext(1,j,k)=zs_ext(1,j,k)
        avgzs_ext(nx,j,k)=zs_ext(nx-1,j,k)
      END DO
    END DO
    DO k=1,nz
      DO j=1,ny
        tem2(1,j,k)=zs(1,j,k)
        tem2(nx,j,k)=zs(nx-1,j,k)
      END DO
    END DO

    IF ( extdopt == 0 ) THEN ! u is staggered as in arps

      CALL avgsu(hgt_ext,nx_ext,ny_ext,nz_ext,1,ny_ext,1,nz_ext,tem5_ext,  &
         tem3_ext)
      tem5_ext(1,1:ny_ext,1:nz_ext) = hgt_ext(1,1:ny_ext,1:nz_ext)
      tem5_ext(nx,1:ny_ext,1:nz_ext) = hgt_ext(nx-1,1:ny_ext,1:nz_ext)
!
      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                     iorder,iprtopt,intropt,iu,ju,xu_ext,yu_ext,            &
                     tem5_ext,avgzs_ext,xu2d,yu2d,tem2,u_ext,               &
                     zsnd,usnd,ubar,u,                                      &
                     dxfldu,dyfldu,rdxfldu,rdyfldu,                         &
                     tem1_ext,tem2_ext,tem3_ext,                            &
                     tem4_ext)

    ELSE  ! u is on a scalar point

      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                     iorder,iprtopt,intropt,iu,ju,x_ext,y_ext,              &
                     hgt_ext,avgzs_ext,xu2d,yu2d,tem2,u_ext,                &
                     zsnd,usnd,ubar,u,                                      &
                     dxfld,dyfld,rdxfld,rdyfld,                             &
                     tem1_ext,tem2_ext,tem3_ext,                            &
                     tem4_ext)

    ENDIF

    DO ksmth=1,nsmooth
      CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                  smfct1,tem2,   u,tem1,   u)
!      call smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,
!  :                smfct2,tem2,   u,tem1,   u)
    END DO

    CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'umin final= ', amin,',  umax final=',amax
!
!-----------------------------------------------------------------------
!
!  Process v component
!
!-----------------------------------------------------------------------
!
    CALL avgsv(zs_ext,nx,ny,nz_ext, 1,nx, 1,nz_ext, avgzs_ext, tem3_ext)
    CALL avgsv(zs,nx,ny,nz, 1,nx, 1,nz, tem2, tem3)

    DO k=1,nz_ext
      DO i=1,nx
        avgzs_ext(i,1,k)=zs_ext(i,1,k)
        avgzs_ext(i,ny,k)=zs_ext(i,ny-1,k)
      END DO
    END DO
    DO k=1,nz
      DO i=1,nx
        tem2(i,1,k)=zs(i,1,k)
        tem2(i,ny,k)=zs(i,ny-1,k)
      END DO
    END DO

    IF ( extdopt == 0 ) THEN ! v is staggered as in arps

      CALL avgsv(hgt_ext,nx_ext,ny_ext,nz_ext,1,nx_ext,1,nz_ext,tem5_ext,  &
         tem3_ext)
      tem5_ext(1:nx_ext,1,1:nz_ext) = hgt_ext(1:nx_ext,1,1:nz_ext)
      tem5_ext(1:nx_ext,ny,1:nz_ext) = hgt_ext(1:nx_ext,ny-1,1:nz_ext)

      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                     iorder,iprtopt,intropt,iv,jv,xv_ext,yv_ext,            &
                     tem5_ext,avgzs_ext,xv2d,yv2d,tem2,v_ext,               &
                     zsnd,vsnd,vbar,v,                                      &
                     dxfldv,dyfldv,rdxfldv,rdyfldv,                         &
                     tem1_ext,tem2_ext,tem3_ext,                            &
                     tem4_ext)

    ELSE  ! v is on a scalar point

      iprtopt=0
      CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                     iorder,iprtopt,intropt,iv,jv,x_ext,y_ext,              &
                     hgt_ext,avgzs_ext,xv2d,yv2d,tem2,v_ext,                &
                     zsnd,vsnd,vbar,v,                                      &
                     dxfld,dyfld,rdxfld,rdyfld,                             &
                     tem1_ext,tem2_ext,tem3_ext,                            &
                     tem4_ext)

    ENDIF

    DO ksmth=1,nsmooth
      CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                  smfct1,tem2,   v,tem1,   v)
    END DO
    CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'vmin final= ', amin,',  vmax final=',amax
!
!-----------------------------------------------------------------------
!
!  Process 2D surface fields if required (sfcout = 1)
!
!-----------------------------------------------------------------------
!
    ! determine the number of soil types supplied
    IF (nstyp_ext <= 0) THEN
      DO is=1,nstyps
        IF (tsfc_ext(1,1,is) > 0) nstyp_ext = is
      ENDDO
    ENDIF

    ! GMB 2001-03-13 wdt
    ! This should never be the case (or else ext sfc arrays allocated to small!):
    nstyp_ext = min(nstyp_ext,nstyps)

    tsfcout = 0
    tsoilout = 0
    wetsout = 0
    wetdout = 0
    wetcout = 0

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!wdt begin block:
!commented out mkarps2d & smooth9's, add call to intrp_soil

    DO is=0,nstyp_ext
      CALL a3dmax0(tsfc_ext(1,1,is),1,nx_ext,1,nx_ext,1,ny_ext,1,ny_ext,    &
           1,1,1,1,amax,amin)
      WRITE(6,'(1x,a,i3,2(a,e13.6))') 'is=',is,                             &
              'tsfc_ext_min= ', amin,', tsfc_ext_max=',amax

      IF ( tsfc_ext(1,1,is) > 0.0 ) THEN
        tsfcout = 1
!      CALL mkarps2d (nx_ext,ny_ext,nx,ny,                                 &
!                     iorder,iscl,jscl,x_ext,y_ext,                        &
!                     xs2d,ys2d,tsfc_ext(1,1,is),tsfc(1,1,is),             &
!                     dxfld,dyfld,rdxfld,rdyfld,                           &
!                     tem1_ext(1,1,1),tem1_ext(1,1,2),                     &
!                     tem1_ext(1,1,3),tem1_ext(1,1,4))
!      DO ksmth=1,nsmooth
!        CALL smooth9p(tsfc(1,1,is), nx,ny, 1,nx,1,ny, tem1)
!      END DO
      END IF

      IF ( tsoil_ext(1,1,is) > 0.0 ) THEN
        tsoilout = 1
!      CALL mkarps2d (nx_ext,ny_ext,nx,ny,                                 &
!                     iorder,iscl,jscl,x_ext,y_ext,                        &
!                     xs2d,ys2d,tsoil_ext(1,1,is),tsoil(1,1,is),           &
!                     dxfld,dyfld,rdxfld,rdyfld,                           &
!                     tem1_ext(1,1,1),tem1_ext(1,1,2),                     &
!                     tem1_ext(1,1,3),tem1_ext(1,1,4))
!      DO ksmth=1,nsmooth
!        CALL smooth9p(tsoil(1,1,is), nx,ny, 1,nx,1,ny, tem1)
!      END DO
      END IF

      IF ( wetsfc_ext(1,1,is) >= 0.0 ) THEN
        wetsout = 1
!      CALL mkarps2d (nx_ext,ny_ext,nx,ny,                                 &
!                     iorder,iscl,jscl,x_ext,y_ext,                        &
!                     xs2d,ys2d,wetsfc_ext(1,1,is),wetsfc(1,1,is),         &
!                     dxfld,dyfld,rdxfld,rdyfld,                           &
!                     tem1_ext(1,1,1),tem1_ext(1,1,2),                     &
!                     tem1_ext(1,1,3),tem1_ext(1,1,4))
!      DO ksmth=1,nsmooth
!        CALL smooth9p(wetsfc(1,1,is), nx,ny, 1,nx,1,ny, tem1)
!      END DO
      END IF

      IF ( wetdp_ext(1,1,is) >= 0.0 ) THEN
        wetdout = 1
!      CALL mkarps2d (nx_ext,ny_ext,nx,ny,                                 &
!                     iorder,iscl,jscl,x_ext,y_ext,                        &
!                     xs2d,ys2d,wetdp_ext(1,1,is),wetdp(1,1,is),           &
!                     dxfld,dyfld,rdxfld,rdyfld,                           &
!                     tem1_ext(1,1,1),tem1_ext(1,1,2),                     &
!                     tem1_ext(1,1,3),tem1_ext(1,1,4))
!      DO ksmth=1,nsmooth
!        CALL smooth9p(wetdp(1,1,is), nx,ny, 1,nx,1,ny, tem1)
!      END DO
      END IF

      IF ( wetcanp_ext(1,1,is) >= 0.0 ) THEN
        wetcout = 1
!      CALL mkarps2d (nx_ext,ny_ext,nx,ny,                                 &
!                     iorder,iscl,jscl,x_ext,y_ext,                        &
!                     xs2d,ys2d,wetcanp_ext(1,1,is),wetcanp(1,1,is),       &
!                     dxfld,dyfld,rdxfld,rdyfld,                           &
!                     tem1_ext(1,1,1),tem1_ext(1,1,2),                     &
!                     tem1_ext(1,1,3),tem1_ext(1,1,4))
!      DO ksmth=1,nsmooth
!        CALL smooth9p(wetcanp(1,1,is), nx,ny, 1,nx,1,ny, tem1)
!      END DO
      END IF
    END DO ! is
!wdt end

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

    IF ( snowdpth_ext(1,1) >= 0 ) THEN

      snowdout = 1

!wdt
!     CALL mkarps2d (nx_ext,ny_ext,nx,ny,              &
!                    iorder,iscl,jscl,x_ext,y_ext,     &
!                    xs2d,ys2d,snowdpth_ext,snowdpth,  &
!                    dxfld,dyfld,rdxfld,rdyfld,        &
!                    tem1_ext(1,1,1),tem1_ext(1,1,2),  &
!                    tem1_ext(1,1,3),tem1_ext(1,1,4))

      DO j=1,ny
        DO i=1,nx
          dmin=((xs2d(i,j)-x_ext(1))*(xs2d(i,j)-x_ext(1)) + &
                (ys2d(i,j)-y_ext(1))*(ys2d(i,j)-y_ext(1)))
          isnow=1
          jsnow=1
          DO jj=jscl(i,j),MAX(jscl(i,j)+1,ny_ext)
            DO ii=iscl(i,j),MAX(iscl(i,j)+1,nx_ext)
              dd=((xs2d(i,j)-x_ext(ii))*(xs2d(i,j)-x_ext(ii)) + &
                  (ys2d(i,j)-y_ext(jj))*(ys2d(i,j)-y_ext(jj)))
              IF(dd < dmin) THEN
                dmin=dd
                isnow=ii
                jsnow=jj
              END IF
            END DO
          END DO
          snowdpth(i,j) = snowdpth_ext(isnow,jsnow)
        END DO
      END DO

    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!wdt begin block:
    xw = 0.
    yw = 0.
    DO j = 1,ny-1
      DO i = 1,nx-1
        xw(i,j) =  &
           MAX(0.,MIN(1.,(x_ext(iscl(i,j)+1)-xs2d(i,j))  &
                      /(x_ext(iscl(i,j)+1)-x_ext(iscl(i,j)))))
        yw(i,j) =  &
           MAX(0.,MIN(1.,(y_ext(jscl(i,j)+1)-ys2d(i,j))  &
                      /(y_ext(jscl(i,j)+1)-y_ext(jscl(i,j)))))
      END DO
    END DO
    !tmp wdt - for reference
    !CALL intrp_soil(nx,ny,nx1,ny1,nstyp,nstyp,xw,yw,ix,jy, &
    !                tsfc,tsoil,wetsfc,wetdp,wetcanp,  &
    !                soiltyp,stypfrct,vegtyp, &
    !                tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,  &
    !                soiltyp1,stypfrct1,vegtyp1)
    CALL intrp_soil(nx_ext,ny_ext,nx,ny,nstyp_ext,nstyp,xw,yw,iscl,jscl, &
                    tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,  &
                    soiltyp_ext,stypfrct_ext,vegtyp_ext, &
                    tsfc,tsoil,wetsfc,wetdp,wetcanp,  &
                    soiltyp,stypfrct,vegtyp)

!wdt FIXME other data types ...; make soiltyp consistant with 
! soil temps, etc.; move this section before  soil temp/moist section
!  INTEGER, allocatable :: soiltyp (:,:,:)! Soil type
!  REAL, allocatable ::    stypfrct(:,:,:)! Soil type fraction
!  INTEGER, allocatable :: vegtyp (:,:)   ! Vegetation type
!  REAL, allocatable ::    lai    (:,:)   ! Leaf Area Index
!  REAL, allocatable ::    roufns (:,:)   ! Surface roughness
!  REAL, allocatable ::    veg    (:,:)   ! Vegetation fraction
!
!  CALL fix_stypfrct_nstyp(nx,ny,nstyp_ext,nstyp,stypfrct)
    
!wdt FIXME
! Add capability to write out sfcdata file here? ...
!wdt end

    IF ( tsfcout == 1 .OR. tsoilout == 1 .OR.                             &
           wetsout == 1 .OR. wetdout == 1 .OR.                            &
           wetcout == 1 .OR. snowdout == 1 ) THEN

      CALL cvttsnd( curtim, timsnd, tmstrln )

      soiloutfl = runname(1:lfnkey)//".soilvar."//timsnd(1:tmstrln)
      lfn = lfnkey + 9 + tmstrln

      IF( dirname /= ' ' ) THEN

        temchar = soiloutfl
        soiloutfl = dirname(1:ldirnam)//'/'//temchar
        lfn  = lfn + ldirnam + 1

      END IF

      IF (soildmp > 0) THEN

        CALL fnversn(soiloutfl, lfn)

        PRINT *, 'Write soil initial data to ',soiloutfl(1:lfn)

        CALL wrtsoil(nx,ny,nstyp, soiloutfl,dx,dy,                       &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,       &
                   tsfcout,tsoilout,wetsout,                              &
                   wetdout,wetcout,snowdout,                              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)

        IF (soildmp == 1) CALL soilcntl(nx,ny, soiloutfl,                 &
                    tsfcout,tsoilout,wetsout,                             &
                    wetdout,wetcout,snowdout,                             &
                    x,y)

      END IF

    END IF
!
!-----------------------------------------------------------------------
!
!  Test code, for diagnostic testing.
!  Find x,y of diagnostic sounding location in ARPS grid.
!
!-----------------------------------------------------------------------
!
    CALL lltoxy(1,1,latdiag,londiag,xdiag,ydiag)
!
    dmin=((xdiag-xscl(1))*(xdiag-xscl(1))+                                &
          (ydiag-yscl(1))*(ydiag-yscl(1)))
    idiag=1
    jdiag=1
!
    DO j=2,ny-2
      DO i=2,nx-2
        dd=((xdiag-xscl(i))*(xdiag-xscl(i))+                              &
            (ydiag-yscl(j))*(ydiag-yscl(j)))
        IF(dd < dmin) THEN
          dmin=dd
          idiag=i
          jdiag=j
        END IF
      END DO
    END DO
    CALL xytoll(1,1,xscl(idiag),yscl(jdiag),                              &
                latd,lond)
    WRITE(6,'(a,f10.4,f10.4,/a,i5,i5,a,f10.4,f10.4)')                     &
          ' Nearest ARPS pt to diagnostic lat,lon: ',                     &
          latdiag,londiag,                                                &
          ' Diagnostic i,j: ',                                            &
          idiag,jdiag,' lat,lon= ',latd,lond
    WRITE(6,'(///a,/2x,a)')                                               &
        ' ARPS extracted sounding at idiag,jdiag',                        &
        'k   pres   hgt   temp   theta   dewp     u     v     dir    spd'
!
!-----------------------------------------------------------------------
!
!  Convert units of ARPS data and write as a sounding.
!
!-----------------------------------------------------------------------
!
    DO k=nz-2,1,-1
      ppasc=pbar(idiag,jdiag,k)+pprt(idiag,jdiag,k)
      pmb=.01*(pbar(idiag,jdiag,k)+pprt(idiag,jdiag,k))
      theta=ptbar(idiag,jdiag,k)+ptprt(idiag,jdiag,k)
      tc=(theta*((ppasc/p0)**rddcp))-273.15
      IF( qv(idiag,jdiag,k) > 0.) THEN
        smix=qv(idiag,jdiag,k)/(1.-qv(idiag,jdiag,k))
        e=(pmb*smix)/(0.62197 + smix)
        bige=e/( 1.001 + ( (pmb - 100.) / 900.) * 0.0034)
        alge = ALOG(bige/6.112)
        tdc = (alge * 243.5) / (17.67 - alge)
      ELSE
        tdc = tc-30.
      END IF

      CALL uvrotdd(1,1,londiag,                                           &
                   u(idiag,jdiag,k),                                      &
                   v(idiag,jdiag,k),                                      &
                   dir,spd)

      WRITE(6,'(i4,f6.0,f9.0,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1,f7.1)')        &
              k,pmb,                                                      &
              zs(idiag,jdiag,k),                                          &
              tc,theta,tdc,                                               &
              u(idiag,jdiag,k),                                           &
              v(idiag,jdiag,k),                                           &
              dir,spd
    END DO
!
!-----------------------------------------------------------------------
!
!  Find vertical velocity and make any u-v adjustments
!  according to wndadj option.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          rhostr(i,j,k)=ABS(j3(i,j,k))*rhobar(i,j,k)
        END DO
      END DO
    END DO

    IF(wndadj == 0) THEN

      IF ( extdopt == 0 ) THEN  ! ARPS history data

        iprtopt=0
        CALL mkarpsvar(nx_ext,ny_ext,nz_ext,nx,ny,nz,lvlprof,                 &
                       iorder,iprtopt,intropt,iscl,jscl,x_ext,y_ext,          &
                       zp_ext,zp2_ext,xs2d,ys2d,zp,w_ext,                     &
                       zsnd,dumsnd,wbar,w,                                    &
                       dxfld,dyfld,rdxfld,rdyfld,                             &
                       tem1_ext,tem2_ext,tem3_ext,                            &
                       tem4_ext)

        DO ksmth=1,nsmooth
          CALL smooth3d(nx,ny,nz, 1,nx,1,ny,1,nz,                             &
                        smfct1,zp,w,tem1,w)
        END DO

      ELSE

        w = 0.

      ENDIF

    ELSE

      CALL adjuvw( nx,ny,nz, u,v,w,wcont,ptprt,ptbar,                     &
                   zp,j1,j2,j3,aj3z,mapfct,rhostr,tem1,                   &
                   wndadj,obropt,obrzero,0,                           &
                   tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

    END IF
!
!-----------------------------------------------------------------------
!
!  Enforce vertical w boundary conditions
!
!-----------------------------------------------------------------------
!
    IF (ext_lbc == 1 .or. ext_vbc == 1)                                   &
       CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)
!
    IF (ext_vbc == 1) CALL vbcw(nx,ny,nz,w,wcont,tbc,bbc,u,v,  &
              rhostr,tem1,tem2,tem3,j1,j2,j3)
!
!
!-----------------------------------------------------------------------
!
!  Assign zero-gradient horizontal boundary conditions
!  to the horizontal & vertical winds.
!
!-----------------------------------------------------------------------
!
    IF (ext_lbc == 1) THEN
      CALL lbcw(nx,ny,nz,1.0, w,wcont,tem1,tem2,tem3,tem4,3,3,3,3)
      CALL bcu(nx,ny,nz,1.0, u, tem1,tem2,tem3,tem4, 3,3,3,3,1,1)
      CALL bcv(nx,ny,nz,1.0, v, tem1,tem2,tem3,tem4, 3,3,3,3,1,1)
    ENDIF
!
!-----------------------------------------------------------------------
!
!  Zero out uninitialized fields
!
!-----------------------------------------------------------------------
!
     tem1=0.
!
!-----------------------------------------------------------------------
!
!  Print out the max/min of output variables.
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/1x,a/)')                                                   &
        'Min and max of External data interpolated to the ARPS grid:'

    CALL a3dmax0(x,1,nx,1,nx,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(/1x,2(a,e13.6))')                                           &
            'xmin    = ', amin,',  xmax    =',amax

    CALL a3dmax0(y,1,ny,1,ny,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'ymin    = ', amin,',  ymax    =',amax

    CALL a3dmax0(z,1,nz,1,nz,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'zmin    = ', amin,',  zmax    =',amax

    CALL a3dmax0(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,                    &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'zpmin   = ', amin,', zpmax    =',amax

    CALL a3dmax0(rhobar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,              &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'rhobarmin=', amin,', rhobarmax=',amax

    CALL a3dmax0(u,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,                     &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'umin    = ', amin,',  umax    =',amax

    CALL a3dmax0(v,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,                     &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'vmin    = ', amin,',  vmax    =',amax

    CALL a3dmax0(w,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,                     &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'wmin    = ', amin,',  wmax    =',amax

    CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,               &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'ptbarmin= ', amin,',  ptbarmax=',amax

    CALL a3dmax0(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,               &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'ptprtmin= ', amin,',  ptprtmax=',amax

    CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,                &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'pbarmin= ', amin,',  pbarmax =',amax

    CALL a3dmax0(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,                &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'pprtmin = ', amin,',  pprtmax =',amax

    CALL a3dmax0(qv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,                  &
              amax,amin)
    WRITE(6,'(1x,2(a,e13.6))')                                            &
            'qvmin   = ', amin,',  qvmax   =',amax

    IF(sfcout == 1) THEN

      CALL a3dmax0(tsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')                                          &
              'tsfcmin = ', amin,',  tsfcmax =',amax
      CALL a3dmax0(tsoil,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')                                          &
              'tsoilmin= ', amin,',  tsoilmax=',amax
      CALL a3dmax0(wetsfc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')                                          &
              'wetsfcmin = ', amin,',  wetsfcmax =',amax
      CALL a3dmax0(wetdp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')                                          &
              'wetdpmin= ', amin,',  wetdpmax=',amax
      CALL a3dmax0(wetcanp,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')                                          &
              'wetcanpmin = ', amin,',  wetcanpmax =',amax

    END IF
!
!-----------------------------------------------------------------------
!
!  Print the mean sounding that was used in setting the
!  mean ARPS variables.
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(a/a)')                                                      &
        '  Mean sounding for ARPS base state variables',                  &
        '  k     p(mb)     z(m)    pt(mb)    T(C)   qv(kg/kg) ',          &
        '  RH %  u(m/s)    v(m/s)'

    DO k=lvlprof,1,-50
      pres = psnd(k)
      temp = ptsnd(k)*((pres/p0)**rddcp)
      rh=AMAX1(0.01,rhmax-(rhssnd(k)*rhssnd(k)))
      qvsat=f_qvsat( pres, temp )
      qvval=rh*qvsat
      WRITE(6,'(i4,f9.1,f9.1,f9.1,f9.1,f9.5,f9.1,f9.1,f9.1)')             &
          k,(0.01*psnd(k)),zsnd(k),ptsnd(k),(temp-273.15),                &
          qvval,(100.*rh),usnd(k),vsnd(k)
    END DO
!
!-----------------------------------------------------------------------
!
!  Data dump of the model grid and base state arrays:
!
!  First find a unique name basdmpfn(1:lbasdmpf) for the grid and
!  base state array dump file
!
!-----------------------------------------------------------------------
!
    tem1=0.

    ldirnam=LEN(dirname)
    CALL strlnth( dirname, ldirnam)

    IF ( hdmpfmt == 5 ) THEN
      WRITE (6,'(a/a)')                                                   &
          'Program ext2arps does not support Savi3D format.',             &
          'Reset to binary format, hdmpfmt=1.'
      hdmpfmt = 1
    END IF

    IF ( hdmpfmt == 9 ) GO TO 450

    CALL gtbasfn(runname(1:lfnkey),dirname,ldirnam,hdmpfmt,               &
        mgrid,nestgrd,basdmpfn,lbasdmpf)

    PRINT*,'Writing base state history dump ',basdmpfn(1:lbasdmpf)

    grdbas =1
    mstout =1
    rainout=0
    prcout =0
    trbout =0
    IF (nstyp < 2) landout = 0

    CALL dtadump(nx,ny,nz,nstyp,                                         &
             hdmpfmt,iniotfu,basdmpfn(1:lbasdmpf),grdbas,filcmprs,        &
             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,                          &
             tem1,tem1,tem1,                                              &
             ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,                      &
             x,y,z,zp,hterain,j1,j2,j3,                                   &
             soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
             tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
             tem1,tem1,tem1,                                              &
             tem1,tem1,tem1,                                              &
             tem1,tem1,tem1,tem1,                                         &
             tem2,tem3,tem4)

    450 CONTINUE

    PRINT*,'curtim = ', curtim

    CALL gtdmpfn(runname(1:lfnkey),dirname,ldirnam,curtim,hdmpfmt,        &
               mgrid,nestgrd, hdmpfn, ldmpf)


    WRITE(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)

    grdbas=0
    CALL dtadump(nx,ny,nz,nstyp,                                         &
             hdmpfmt,iniotfu,hdmpfn(1:ldmpf),grdbas,filcmprs,             &
             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,                          &
             tem1,tem1,tem1,                                              &
             ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar,                      &
             x,y,z,zp,hterain,j1,j2,j3,                                   &
             soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
             tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
             tem1,tem1,tem1,                                              &
             tem1,tem1,tem1,                                              &
             tem1,tem1,tem1,tem1,                                         &
             tem2,tem3,tem4)

  END DO
!
!-----------------------------------------------------------------------
!
!  Friendly exit message.
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a/a,i4,a)')                                                 &
          ' Normal succesful completion of EXT2ARPS.',                  &
          ' Processed',nextdfil,' file(s)'
  STOP
!
!-----------------------------------------------------------------------
!
!  Problem doing time conversions.
!
!-----------------------------------------------------------------------
!
  920 CONTINUE
  WRITE(6,'(a,/,a,i4,a,i4/,a,a19)')                                     &
          ' Aborting, error in time format for external file',          &
          ' File number:',ifile,' of',nextdfil,                         &
          ' External file time provided:',extdtime(ifile)
  STOP
!
!-----------------------------------------------------------------------
!
!  Error status returned from rdextfil
!
!-----------------------------------------------------------------------
!
  999 CONTINUE
  WRITE(6,'(a,i6)')                                                     &
          ' Aborting, error reading external file. istatus=',           &
            istatus
  STOP
END PROGRAM ext2arps