PROGRAM arpsdiff,36
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM ARPSDIFF                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Computes difference between two ARPS history files and writes
!  the difference as a similar history file.
!
!  Reads in a history file produced by ARPS in any ARPS format.
!
!  Parameters grdout,varout,mstout,iceout and trbout should be input
!  with the same values as in the data dump subroutines in the model.
!
!-----------------------------------------------------------------------
!
!
!  AUTHOR: Keith Brewster  OU School of Meteorology. February 1992
!
!  MODIFICATION HISTORY:
!   11 Aug 1992  (KB) changed from arps2.5 to arps3.0
!   19 May 1992  (KB) changed from arps3.0 to arps3.2
!   24 May 1992  (KB) added code to handle staggered grid
!                     e.g., intrpu,intrpv,intrpw calls
!   12 Oct 1993  (KB) modified storage to save memory space
!                Differences are now output to interpolated arrays.
!
!   08/01/1995  (Ming Xue)
!   Changed the difference to be that between two total fields.
!   The output perturbation fields represent the difference in
!   the total fields.  The mean-state fields are not changed.
!
!   09/07/1995  (KB)
!   Added differencing and writing of surface (soil) fields.
!   This change affects the input file.
!
!   09/12/1995 (KB)
!   Cleared some bugs and restored some diagnostic printing.
!
!   03/13/1996 (Ming Xue)
!   Added array tke, atke, vtke. Change km to kmh and kmv.
!
!   04/02/1996 (Keith Brewster)
!   New runname is now read-in instead of output file names.
!   File names are constructed from runname just as in arps.
!
!   04/30/1996 (Keith Brewster)
!   Upgraded interpolation to n-order (1-to-fourth) polynomial
!   based on Gauss forward method.  Required some reorganization
!   to use the routine efficiently.  Added "iorder" as input
!   variable.
!
!   11/07/1996 (Keith Brewster)
!   Modified code for new interpolation scheme.
!
!   12/16/1996 (Yuhe Liu)
!   Corrected the dimension definitions for all verification of soil
!   and vegetation variables. They should be (vnx,vny) instead of
!   (nx,ny).
!
!   04/23/1998 (Yvette Richardson)
!   Corrected a bug affecting results when fdx.ne.fdy when
!   calculating fy0.
!
!   12/14/1998 (Donghai Wang)
!   Added the snow cover.
!
!   10/15/1999 (KB via Eric Kemp)
!   Corrected dimension statement of vxs and vys.  Affects only
!   calculation of vxs(vnx) and vys(vny).
!
!   11/03/1999 (KB via Eric Kemp)
!   Corrected dimension of snow cover and sfc flux variables
!   that had recently been added.
!
!   7 May 2002 (Eric Kemp)
!   Added allocation for array lon, and set nstyp = nstyps to 
!   correctly read in soil data.  Also, added code to skip interpolation 
!   when both grids are identical.
!
!-----------------------------------------------------------------------
!
!  DATA ARRAYS READ IN:
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    u        x component of velocity (m/s)
!    v        y component of velocity (m/s)
!    w        vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qv       water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  CALCULATED DATA ARRAYS:
!
!    uprt    perturbation x component of velocity (m/s)
!    vprt    perturbation y component of velocity (m/s)
!    wprt    perturbation z component of velocity (m/s)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx       ! Number of grid points in the x-direction
  INTEGER :: ny       ! Number of grid points in the y-direction
  INTEGER :: nz       ! Number of grid points in the z-direction
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
!  Arrays to be read in:
!
!-----------------------------------------------------------------------
!
  INTEGER :: nstyps            ! Maximum number of soil types in each
                               ! grid box
  PARAMETER (nstyps=4)

  REAL, ALLOCATABLE :: x    (:)          ! The x-coord. of the physical and
                                         ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: y    (:)          ! The y-coord. of the physical and
                                         ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: z    (:)          ! The z-coord. of the computational grid.
                                         ! Defined at w-point on the staggered grid.
  REAL, ALLOCATABLE :: zp   (:,:,:)      ! The physical height coordinate defined at
                                         ! w-point of the staggered grid.

  REAL, ALLOCATABLE :: j1     (:,:,:)    ! Coordinate transformation Jacobian defined
                                         ! as - d( zp )/d( x )
  REAL, ALLOCATABLE :: j2     (:,:,:)    ! Coordinate transformation Jacobian defined
                                         ! as - d( zp )/d( y )
  REAL, ALLOCATABLE :: j3     (:,:,:)    ! Coordinate transformation Jacobian defined
                                         ! as d( zp )/d( z )
  REAL, ALLOCATABLE :: hterain(:,:)      ! Terrain height.

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

  REAL, ALLOCATABLE :: tke    (:,:,:)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: kmh    (:,:,:)    ! Horizontal turb. mixing coef. for
                                         ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: kmv    (:,:,:)    ! Vertical turb. mixing coef. for
                                         ! momentum. ( m**2/s )

  REAL, ALLOCATABLE :: ubar   (:,:,:)    ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vbar   (:,:,:)    ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: wbar   (:,:,:)    ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: ptbar  (:,:,:)    ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: pbar   (:,:,:)    ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: rhobar (:,:,:)    ! Base state air density (kg/m**3)
  REAL, ALLOCATABLE :: qvbar  (:,:,:)    ! Base state water vapor specific humidity

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

  REAL, ALLOCATABLE :: tsfc   (:,:,:)    ! Temperature at surface (K)
  REAL, ALLOCATABLE :: tsoil  (:,:,:)    ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: wetsfc (:,:,:)    ! Surface soil moisture
  REAL, ALLOCATABLE :: wetdp  (:,:,:)    ! Deep soil moisture
  REAL, ALLOCATABLE :: wetcanp(:,:,:)    ! Canopy water amount
  REAL, ALLOCATABLE :: snowdpth(:,:)     ! Snow depth (m)

  REAL, ALLOCATABLE :: raing  (:,:)      ! Cumulus convective rain
  REAL, ALLOCATABLE :: rainc  (:,:)      ! Cumulus convective rain
  REAL, ALLOCATABLE :: prcrate(:,:,:)    ! precipitation rate (kg/(m**2*s))
                                         ! prcrate(1,1,1) = total precip. rate
                                         ! prcrate(1,1,2) = grid scale precip. rate
                                         ! prcrate(1,1,3) = cumulus precip. rate
                                         ! prcrate(1,1,4) = microphysics precip. rate

  REAL, ALLOCATABLE :: radfrc(:,:,:)     ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: radsw (:,:)       ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: rnflx (:,:)       ! Net radiation flux absorbed by surface

  REAL, ALLOCATABLE :: usflx (:,:)       ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: vsflx (:,:)       ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: ptsflx(:,:)       ! Surface heat flux (K*kg/(m*s**2))
  REAL, ALLOCATABLE :: qvsflx(:,:)       ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Verification Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: vx     (:)     ! The x-coord. of the physical and
                                      ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: vy     (:)     ! The y-coord. of the physical and
                                      ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: vz     (:)     ! The z-coord. of the computational grid.
                                      ! Defined at w-point on the staggered grid.
  REAL, ALLOCATABLE :: vzp    (:,:,:) ! The physical height coordinate defined at
                                      ! w-point of the staggered grid.
  REAL, ALLOCATABLE :: vuprt  (:,:,:) ! Perturbation u-velocity (m/s)
  REAL, ALLOCATABLE :: vvprt  (:,:,:) ! Perturbation v-velocity (m/s)
  REAL, ALLOCATABLE :: vwprt  (:,:,:) ! Perturbation w-velocity (m/s)
  REAL, ALLOCATABLE :: vptprt (:,:,:) ! Perturbation potential temperature (K)
  REAL, ALLOCATABLE :: vpprt  (:,:,:) ! Perturbation pressure (Pascal)
  REAL, ALLOCATABLE :: vqvprt (:,:,:) ! Perturbation water vapor specific humidity
  REAL, ALLOCATABLE :: vqc    (:,:,:) ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqr    (:,:,:) ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqi    (:,:,:) ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqs    (:,:,:) ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqh    (:,:,:) ! Hail mixing ratio (kg/kg)

  REAL, ALLOCATABLE :: vtke   (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: vkmh   (:,:,:) ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: vkmv   (:,:,:) ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: vubar  (:,:,:) ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vvbar  (:,:,:) ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: vwbar  (:,:,:) ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: vptbar (:,:,:) ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: vpbar  (:,:,:) ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: vrhobar(:,:,:) ! Base state air density (kg/m**3)
  REAL, ALLOCATABLE :: vqvbar (:,:,:) ! Base state water vapor specific humidity

  INTEGER, ALLOCATABLE :: vsoiltyp (:,:,:)   ! Soil type
  REAL, ALLOCATABLE :: vstypfrct(:,:,:)      ! Soil type
  INTEGER, ALLOCATABLE :: vvegtyp  (:,:)     ! Vegetation type
  REAL, ALLOCATABLE :: vlai     (:,:)        ! Leaf Area Index
  REAL, ALLOCATABLE :: vroufns  (:,:)        ! Surface roughness
  REAL, ALLOCATABLE :: vveg     (:,:)        ! Vegetation fraction

  REAL, ALLOCATABLE :: vtsfc   (:,:,:)     ! Temperature at surface (K)
  REAL, ALLOCATABLE :: vtsoil  (:,:,:)     ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: vwetsfc (:,:,:)     ! Surface soil moisture
  REAL, ALLOCATABLE :: vwetdp  (:,:,:)     ! Deep soil moisture
  REAL, ALLOCATABLE :: vwetcanp(:,:,:)     ! Canopy water amount
  REAL, ALLOCATABLE :: vsnowdpth(:,:)      ! Snow depth (m)

  REAL, ALLOCATABLE :: vraing(:,:)         ! Grid supersaturation rain
  REAL, ALLOCATABLE :: vrainc(:,:)         ! Cumulus convective rain
  REAL, ALLOCATABLE :: vprcrate(:,:,:)     ! precipitation rate (kg/(m**2*s))
                                           ! prcrate(1,1,1) = total precip. rate
                                           ! prcrate(1,1,2) = grid scale precip. rate
                                           ! prcrate(1,1,3) = cumulus precip. rate
                                           ! prcrate(1,1,4) = microphysics precip. rate

  REAL, ALLOCATABLE :: vradfrc(:,:,:)      ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: vradsw (:,:)        ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: vrnflx (:,:)        ! Net radiation flux absorbed by surface

  REAL, ALLOCATABLE :: vusflx (:,:)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: vvsflx (:,:)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: vptsflx(:,:)        ! Surface heat flux (K*kg/(m*s**2))
  REAL, ALLOCATABLE :: vqvsflx(:,:)        ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Verification data interpolated to model grid
!  These arrays also hold difference fields after call to diffgr.
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: auprt  (:,:,:)    ! Perturbation u-velocity (m/s)
  REAL, ALLOCATABLE :: avprt  (:,:,:)    ! Perturbation v-velocity (m/s)
  REAL, ALLOCATABLE :: awprt  (:,:,:)    ! Perturbation w-velocity (m/s)
  REAL, ALLOCATABLE :: aptprt (:,:,:)    ! Perturbation potential temperature (K)
  REAL, ALLOCATABLE :: apprt  (:,:,:)    ! Perturbation pressure (Pascal)
  REAL, ALLOCATABLE :: aqvprt (:,:,:)    ! Perturbation water vapor specific humidity
  REAL, ALLOCATABLE :: aqc    (:,:,:)    ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqr    (:,:,:)    ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqi    (:,:,:)    ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqs    (:,:,:)    ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: aqh    (:,:,:)    ! Hail mixing ratio (kg/kg)

  REAL, ALLOCATABLE :: atke   (:,:,:)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: akmh   (:,:,:)    ! Horizontal turb. mixing coef. for
                                         ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: akmv   (:,:,:)    ! Vertical turb. mixing coef. for
                                         ! momentum. ( m**2/s )

  REAL, ALLOCATABLE :: aubar  (:,:,:)    ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: avbar  (:,:,:)    ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: awbar  (:,:,:)    ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: aptbar (:,:,:)    ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: apbar  (:,:,:)    ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: arhobar(:,:,:)    ! Base state air density (kg/m**3)
  REAL, ALLOCATABLE :: aqvbar (:,:,:)    ! Base state water vapor specific humidity

  INTEGER, ALLOCATABLE :: asoiltyp (:,:,:)   ! Soil type
  REAL, ALLOCATABLE :: astypfrct(:,:,:)      ! Soil type
  INTEGER, ALLOCATABLE :: avegtyp  (:,:)     ! Vegetation type
  REAL, ALLOCATABLE :: alai     (:,:)        ! Leaf Area Index
  REAL, ALLOCATABLE :: aroufns  (:,:)        ! Surface roughness
  REAL, ALLOCATABLE :: aveg     (:,:)        ! Vegetation fraction

  REAL, ALLOCATABLE :: atsfc   (:,:,:)     ! Temperature at surface (K)
  REAL, ALLOCATABLE :: atsoil  (:,:,:)     ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: awetsfc (:,:,:)     ! Surface soil moisture
  REAL, ALLOCATABLE :: awetdp  (:,:,:)     ! Deep soil moisture
  REAL, ALLOCATABLE :: awetcanp(:,:,:)     ! Canopy water amount
  REAL, ALLOCATABLE :: asnowdpth(:,:)      ! Snow depth (m)

  REAL, ALLOCATABLE :: araing (:,:)        ! Grid supersaturation rain
  REAL, ALLOCATABLE :: arainc (:,:)        ! Cumulus convective rain
  REAL, ALLOCATABLE :: aprcrate(:,:,:)     ! precipitation rate (kg/(m**2*s))
                                           ! prcrate(1,1,1) = total precip. rate
                                           ! prcrate(1,1,2) = grid scale precip. rate
                                           ! prcrate(1,1,3) = cumulus precip. rate
                                           ! prcrate(1,1,4) = microphysics precip. rate

  REAL, ALLOCATABLE :: aradfrc(:,:,:)     ! Radiation forcing (K/s)
  REAL, ALLOCATABLE :: aradsw (:,:)       ! Solar radiation reaching the surface
  REAL, ALLOCATABLE :: arnflx (:,:)       ! Net radiation flux absorbed by surface

  REAL, ALLOCATABLE :: ausflx (:,:)       ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: avsflx (:,:)       ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, ALLOCATABLE :: aptsflx(:,:)       ! Surface heat flux (K*kg/(m*s**2))
  REAL, ALLOCATABLE :: aqvsflx(:,:)       ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)

  REAL, ALLOCATABLE :: vtem1(:,:,:)
  REAL, ALLOCATABLE :: vtem2(:,:,:)
  REAL, ALLOCATABLE :: vtem3(:,:,:)

  REAL, ALLOCATABLE :: xs(:)
  REAL, ALLOCATABLE :: ys(:)
  REAL, ALLOCATABLE :: zps(:,:,:)
  REAL, ALLOCATABLE :: x2d(:,:)
  REAL, ALLOCATABLE :: y2d(:,:)
  REAL, ALLOCATABLE :: lat(:,:),lon(:,:)

  REAL, ALLOCATABLE :: vxs(:)
  REAL, ALLOCATABLE :: vys(:)
  REAL, ALLOCATABLE :: vzps(:,:,:)
  REAL, ALLOCATABLE :: dxfld(:)
  REAL, ALLOCATABLE :: dyfld(:)
  REAL, ALLOCATABLE :: rdxfld(:)
  REAL, ALLOCATABLE :: rdyfld(:)

  INTEGER, ALLOCATABLE :: iloc(:,:),jloc(:,:)
  REAL, ALLOCATABLE :: zpver(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: filename,vfilename,grdbasfn,vgrdbasfn,fcrnam,runnmin
  INTEGER :: lengbf,lenfil
  INTEGER :: ifproj,ivproj
  REAL :: flatnot(2),vlatnot(2)
  REAL :: fscale,ftrulon,fdx,fdy,fx0,fy0
  REAL :: fctrlat,fctrlon
  REAL :: vscale,vtrulon,vdx,vdy,vx0,vy0
  REAL :: vctrlat,vctrlon
  REAL :: time,xctr,yctr
  INTEGER :: i,j,k
  INTEGER :: grdbas
  INTEGER :: iorder,hinfmt,vhinfmt
  INTEGER :: ireturn
  LOGICAL :: comcoord,comcoord2
  INTEGER :: nch

!-----------------------------------------------------------------------
!  The following defines the model dimension parameters for the
!  verification data.
!-----------------------------------------------------------------------
!
  INTEGER :: vnx       ! Number of grid points in the x-direction
  INTEGER :: vny       ! Number of grid points in the y-direction
  INTEGER :: vnz       ! Number of grid points in the z-direction


  INTEGER :: vnxy, vnxz, vnyz, vnxyz

!-----------------------------------------------------------------------
  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  namelist Declarations:
!
!-----------------------------------------------------------------------
!
  NAMELIST /grid_dims/ nx,ny,nz

  NAMELIST /vgrid_dims/ vnx,vny,vnz

  NAMELIST /frcst_fn/ iorder, hinfmt, grdbasfn, filename

  NAMELIST /vrftn_fn/ vhinfmt, vgrdbasfn, vfilename

  NAMELIST /output/ runnmin, hdmpfmt,                      &
                    grdout,   basout,    varout,           &
                    mstout,   iceout,    trbout,           &
                    sfcout,   rainout,   prcout,           &
                    radout,   flxout,    filcmprs

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  nstyp = nstyps ! EMK

  WRITE(6,'(/6(/5x,a))')                                                &
      '###############################################################',&
      '#                                                             #',&
      '# Welcome to ARPSDIFF, a program that reads in history files  #',&
      '# generated by ARPS and produces difference grids.            #',&
      '#                                                             #',&
      '###############################################################'
!
!-----------------------------------------------------------------------
!
!  Set up the default values for all the variables to be read in
!  using the namelist method. In case the user does not specify a
!  particular value, this value will be used.
!
!-----------------------------------------------------------------------
!
   nx   = 67
   ny   = 67
   nz   = 35

   vnx   = 163
   vny   = 99
   vnz   = 43

   iorder   = 3
   hinfmt   = 10
   grdbasfn = 'may20.grbgrdbas'
   filename = 'may20.grb003600'

   runnmin = 'diffm20'
   hdmpfmt = 10
   grdout = 0
   basout = 0
   varout = 1
   mstout = 1
   iceout = 0
   trbout = 0
   sfcout = 0
   rainout = 0
   prcout = 0
   radout = 0
   flxout = 0
   filcmprs = 0

!-----------------------------------------------------------------------
!
!  Read in nx,ny,nz and vnx, vny, vnz
!
!----------------------------------------------------------------------

   READ(5,grid_dims,END=100)
   WRITE(6,'(a)')'Namelist block grid_dims sucessfully read in.' 

   READ(5,vgrid_dims,END=100)
   WRITE(6,'(a)')'Namelist block vgrid_dims sucessfully read in.' 

   vnxy=vnx*vny
   vnxz=vnx*vnz
   vnyz=vny*vnz
   vnxyz=vnxy*vnz
!
!-----------------------------------------------------------------------
!  Allocate and Initialize the variables
!----------------------------------------------------------------------
 
  istatus=0

  ALLOCATE(  x=0
  ALLOCATE(  y=0
  ALLOCATE(  z=0
  ALLOCATE(zp(nx,ny,nz),STAT=istatus)
  zp=0
  ALLOCATE(j1(nx,ny,nz),STAT=istatus)
  j1=0
  ALLOCATE(j2(nx,ny,nz),STAT=istatus)
  j2=0
  ALLOCATE(j3(nx,ny,nz),STAT=istatus)
  j3=0
  ALLOCATE(hterain(nx,ny),STAT=istatus)
  hterain=0
  ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
  uprt=0
  ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
  vprt=0
  ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
  wprt=0
  ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
  ptprt=0
  ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
  pprt=0
  ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
  qvprt=0
  ALLOCATE(qc(nx,ny,nz),STAT=istatus)
  qc=0
  ALLOCATE(qr(nx,ny,nz),STAT=istatus)
  qr=0
  ALLOCATE(qi(nx,ny,nz),STAT=istatus)
  qi=0
  ALLOCATE(qs(nx,ny,nz),STAT=istatus)
  qs=0
  ALLOCATE(qh(nx,ny,nz),STAT=istatus)
  qh=0
  ALLOCATE(tke(nx,ny,nz),STAT=istatus)
  tke=0
  ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
  kmh=0
  ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
  kmv=0
  ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
  ubar=0
  ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
  vbar=0
  ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
  wbar=0
  ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
  ptbar=0
  ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
  pbar=0
  ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
  rhobar=0
  ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
  qvbar=0
  ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
  soiltyp=0
  ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
  stypfrct=0
  ALLOCATE(vegtyp(nx,ny),STAT=istatus)
  vegtyp=0
  ALLOCATE(lai(nx,ny),STAT=istatus)
  lai=0
  ALLOCATE(roufns(nx,ny),STAT=istatus)
  roufns=0
  ALLOCATE(veg(nx,ny),STAT=istatus)
  veg=0
  ALLOCATE(tsfc(nx,ny,0:nstyps),STAT=istatus)
  tsfc=0
  ALLOCATE(tsoil(nx,ny,0:nstyps),STAT=istatus)
  tsoil=0
  ALLOCATE(wetsfc(nx,ny,0:nstyps),STAT=istatus)
  wetsfc=0
  ALLOCATE(wetdp(nx,ny,0:nstyps),STAT=istatus)
  wetdp=0
  ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
  wetcanp=0
  ALLOCATE(snowdpth(nx,ny),STAT=istatus)
  snowdpth=0
  ALLOCATE(raing(nx,ny),STAT=istatus)
  raing=0
  ALLOCATE(rainc(nx,ny),STAT=istatus)
  rainc=0
  ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
  prcrate=0
  ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
  radfrc=0
  ALLOCATE(radsw(nx,ny),STAT=istatus)
  radsw=0
  ALLOCATE(rnflx(nx,ny),STAT=istatus)
  rnflx=0
  ALLOCATE(usflx(nx,ny),STAT=istatus)
  usflx=0
  ALLOCATE(vsflx(nx,ny),STAT=istatus)
  vsflx=0
  ALLOCATE(ptsflx(nx,ny),STAT=istatus)
  ptsflx=0
  ALLOCATE(qvsflx(nx,ny),STAT=istatus)
  qvsflx=0
  ALLOCATE(vx(vnx),STAT=istatus)
  vx=0
  ALLOCATE(vy(vny),STAT=istatus)
  vy=0
  ALLOCATE(vz(vnz),STAT=istatus)
  vz=0
  ALLOCATE(vzp(vnx,vny,vnz),STAT=istatus)
  vzp=0
  ALLOCATE(vuprt(vnx,vny,vnz),STAT=istatus)
  vuprt=0
  ALLOCATE(vvprt(vnx,vny,vnz),STAT=istatus)
  vvprt=0
  ALLOCATE(vwprt(vnx,vny,vnz),STAT=istatus)
  vwprt=0
  ALLOCATE(vptprt(vnx,vny,vnz),STAT=istatus)
  vptprt=0
  ALLOCATE(vpprt(vnx,vny,vnz),STAT=istatus)
  vpprt=0
  ALLOCATE(vqvprt(vnx,vny,vnz),STAT=istatus)
  vqvprt=0
  ALLOCATE(vqc(vnx,vny,vnz),STAT=istatus)
  vqc=0
  ALLOCATE(vqr(vnx,vny,vnz),STAT=istatus)
  vqr=0
  ALLOCATE(vqi(vnx,vny,vnz),STAT=istatus)
  vqi=0
  ALLOCATE(vqs(vnx,vny,vnz),STAT=istatus)
  vqs=0
  ALLOCATE(vqh(vnx,vny,vnz),STAT=istatus)
  vqh=0
  ALLOCATE(vtke(nx,ny,nz),STAT=istatus)
  vtke=0
  ALLOCATE(vkmh(nx,ny,nz),STAT=istatus)
  vkmh=0
  ALLOCATE(vkmv(nx,ny,nz),STAT=istatus)
  vkmv=0
  ALLOCATE(vubar(vnx,vny,vnz),STAT=istatus)
  vubar=0
  ALLOCATE(vvbar(vnx,vny,vnz),STAT=istatus)
  vvbar=0
  ALLOCATE(vwbar(vnx,vny,vnz),STAT=istatus)
  vwbar=0
  ALLOCATE(vptbar(vnx,vny,vnz),STAT=istatus)
  vptbar=0
  ALLOCATE(vpbar(vnx,vny,vnz),STAT=istatus)
  vpbar=0
  ALLOCATE(vrhobar(vnx,vny,vnz),STAT=istatus)
  vrhobar=0
  ALLOCATE(vqvbar(vnx,vny,vnz),STAT=istatus)
  vqvbar=0
  ALLOCATE(vsoiltyp(vnx,vny,nstyps),STAT=istatus)
  vsoiltyp=0
  ALLOCATE(vstypfrct(vnx,vny,nstyps),STAT=istatus)
  vstypfrct=0
  ALLOCATE(vvegtyp(vnx,vny),STAT=istatus)
  vvegtyp=0
  ALLOCATE(vlai(vnx,vny),STAT=istatus)
  vlai=0
  ALLOCATE(vroufns(vnx,vny),STAT=istatus)
  vroufns=0
  ALLOCATE(vveg(vnx,vny),STAT=istatus)
  vveg=0
  ALLOCATE(vtsfc(vnx,vny,0:nstyps),STAT=istatus)
  vtsfc=0
  ALLOCATE(vtsoil(vnx,vny,0:nstyps),STAT=istatus)
  vtsoil=0
  ALLOCATE(vwetsfc(vnx,vny,0:nstyps),STAT=istatus)
  vwetsfc=0
  ALLOCATE(vwetdp(vnx,vny,0:nstyps),STAT=istatus)
  vwetdp=0
  ALLOCATE(vwetcanp(vnx,vny,0:nstyps),STAT=istatus)
  vwetcanp=0
  ALLOCATE(vsnowdpth(vnx,vny),STAT=istatus)
  vsnowdpth=0
  ALLOCATE(vraing(vnx,vny),STAT=istatus)
  vraing=0
  ALLOCATE(vrainc(vnx,vny),STAT=istatus)
  vrainc=0
  ALLOCATE(vprcrate(vnx,vny,4),STAT=istatus)
  vprcrate=0
  ALLOCATE(vradfrc(vnx,vny,vnz),STAT=istatus)
  vradfrc=0
  ALLOCATE(vradsw(vnx,vny),STAT=istatus)
  vradsw=0
  ALLOCATE(vrnflx(vnx,vny),STAT=istatus)
  vrnflx=0
  ALLOCATE(vusflx(vnx,vny),STAT=istatus)
  vusflx=0
  ALLOCATE(vvsflx(vnx,vny),STAT=istatus)
  vvsflx=0
  ALLOCATE(vptsflx(vnx,vny),STAT=istatus)
  vptsflx=0
  ALLOCATE(vqvsflx(vnx,vny),STAT=istatus)
  vqvsflx=0
  ALLOCATE(auprt(nx,ny,nz),STAT=istatus)
  auprt=0
  ALLOCATE(avprt(nx,ny,nz),STAT=istatus)
  avprt=0
  ALLOCATE(awprt(nx,ny,nz),STAT=istatus)
  awprt=0
  ALLOCATE(aptprt(nx,ny,nz),STAT=istatus)
  aptprt=0
  ALLOCATE(apprt(nx,ny,nz),STAT=istatus)
  apprt=0
  ALLOCATE(aqvprt(nx,ny,nz),STAT=istatus)
  aqvprt=0
  ALLOCATE(aqc(nx,ny,nz),STAT=istatus)
  aqc=0
  ALLOCATE(aqr(nx,ny,nz),STAT=istatus)
  aqr=0
  ALLOCATE(aqi(nx,ny,nz),STAT=istatus)
  aqi=0
  ALLOCATE(aqs(nx,ny,nz),STAT=istatus)
  aqs=0
  ALLOCATE(aqh(nx,ny,nz),STAT=istatus)
  aqh=0
  ALLOCATE(atke(nx,ny,nz),STAT=istatus)
  atke=0
  ALLOCATE(akmh(nx,ny,nz),STAT=istatus)
  akmh=0
  ALLOCATE(akmv(nx,ny,nz),STAT=istatus)
  akmv=0
  ALLOCATE(aubar(nx,ny,nz),STAT=istatus)
  aubar=0
  ALLOCATE(avbar(nx,ny,nz),STAT=istatus)
  avbar=0
  ALLOCATE(awbar(nx,ny,nz),STAT=istatus)
  awbar=0
  ALLOCATE(aptbar(nx,ny,nz),STAT=istatus)
  aptbar=0
  ALLOCATE(apbar(nx,ny,nz),STAT=istatus)
  apbar=0
  ALLOCATE(arhobar(nx,ny,nz),STAT=istatus)
  arhobar=0
  ALLOCATE(aqvbar(nx,ny,nz),STAT=istatus)
  aqvbar=0
  ALLOCATE(asoiltyp(nx,ny,nstyps),STAT=istatus)
  asoiltyp=0
  ALLOCATE(astypfrct(nx,ny,nstyps),STAT=istatus)
  astypfrct=0
  ALLOCATE(avegtyp(nx,ny),STAT=istatus)
  avegtyp=0
  ALLOCATE(alai(nx,ny),STAT=istatus)
  alai=0
  ALLOCATE(aroufns(nx,ny),STAT=istatus)
  aroufns=0
  ALLOCATE(aveg(nx,ny),STAT=istatus)
  aveg=0
  ALLOCATE(atsfc(nx,ny,0:nstyps),STAT=istatus)
  atsfc=0
  ALLOCATE(atsoil(nx,ny,0:nstyps),STAT=istatus)
  atsoil=0
  ALLOCATE(awetsfc(nx,ny,0:nstyps),STAT=istatus)
  awetsfc=0
  ALLOCATE(awetdp(nx,ny,0:nstyps),STAT=istatus)
  awetdp=0
  ALLOCATE(awetcanp(nx,ny,0:nstyps),STAT=istatus)
  awetcanp=0
  ALLOCATE(asnowdpth(nx,ny),STAT=istatus)
  asnowdpth=0
  ALLOCATE(araing(nx,ny),STAT=istatus)
  araing=0
  ALLOCATE(arainc(nx,ny),STAT=istatus)
  arainc=0
  ALLOCATE(aprcrate(nx,ny,4),STAT=istatus)
  aprcrate=0
  ALLOCATE(aradfrc(nx,ny,nz),STAT=istatus)
  aradfrc=0
  ALLOCATE(aradsw(nx,ny),STAT=istatus)
  aradsw=0
  ALLOCATE(arnflx(nx,ny),STAT=istatus)
  arnflx=0
  ALLOCATE(ausflx(nx,ny),STAT=istatus)
  ausflx=0
  ALLOCATE(avsflx(nx,ny),STAT=istatus)
  avsflx=0
  ALLOCATE(aptsflx(nx,ny),STAT=istatus)
  aptsflx=0
  ALLOCATE(aqvsflx(nx,ny),STAT=istatus)
  aqvsflx=0
  ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
  tem1=0
  ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
  tem2=0
  ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
  tem3=0
  ALLOCATE(vtem1(vnx,vny,vnz),STAT=istatus)
  vtem1=0
  ALLOCATE(vtem2(vnx,vny,vnz),STAT=istatus)
  vtem2=0
  ALLOCATE(vtem3(vnx,vny,vnz),STAT=istatus)
  vtem3=0
  ALLOCATE(xs(nx),STAT=istatus)
  xs=0
  ALLOCATE(ys(ny),STAT=istatus)
  ys=0
  ALLOCATE(zps(nx,ny,nz),STAT=istatus)
  zps=0
  ALLOCATE(x2d(nx,ny),STAT=istatus)
  x2d=0
  ALLOCATE(y2d(nx,ny),STAT=istatus)
  y2d=0
  ALLOCATE(lat(nx,ny),STAT=istatus)
  lat=0
  ALLOCATE(lon(nx,ny),STAT=istatus) ! EMK
  lon=0 ! EMK
  ALLOCATE(vxs(vnx),STAT=istatus)
  vxs=0
  ALLOCATE(vys(vny),STAT=istatus)
  vys=0
  ALLOCATE(vzps(vnx,vny,vnz),STAT=istatus)
  vzps=0
  ALLOCATE(dxfld(vnx),STAT=istatus)
  dxfld=0
  ALLOCATE(dyfld(vny),STAT=istatus)
  dyfld=0
  ALLOCATE(rdxfld(vnx),STAT=istatus)
  rdxfld=0
  ALLOCATE(rdyfld(vny),STAT=istatus)
  rdyfld=0
  ALLOCATE(iloc(nx,ny),STAT=istatus)
  iloc=0
  ALLOCATE(jloc(nx,ny),STAT=istatus)
  jloc=0
  ALLOCATE(zpver(nx,ny,vnz),STAT=istatus)
  zpver=0

  mgrid = 1
  nestgrd = 0
  grbpkbit = 16

  READ(5,frcst_fn,END=100)
  WRITE(6,'(a)')'Namelist block frcst_fn sucessfully read in.' 

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

  lenfil = LEN_trim(filename)
  WRITE(6,'(/a,/1x,a)')' The data set name is ', filename(1:lenfil)
!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
  CALL dtaread(nx,ny,nz,nstyps,                                         &
               hinfmt,nch,grdbasfn(1:lengbf),lengbf,                    &
               filename(1:lenfil),lenfil,time,                          &
               x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,                &
               qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                  &
               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,            &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,                                      &
               usflx,vsflx,ptsflx,qvsflx,                               &
               ireturn, tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
  IF( ireturn == 0 ) THEN   ! successful read
    curtim=time
    fcrnam=runname
    ifproj=mapproj
    fscale=sclfct
    flatnot(1)=trulat1
    flatnot(2)=trulat2
    ftrulon=trulon
    fdx=x(3)-x(2)
    fdy=y(3)-y(2)
    fctrlat=ctrlat
    fctrlon=ctrlon
    CALL setmapr(ifproj,fscale,flatnot,ftrulon)
    CALL lltoxy(1,1,fctrlat,fctrlon,xctr,yctr)
    fx0=xctr-fdx*((nx-3)/2)
    fy0=yctr-fdy*((ny-3)/2)
    CALL setorig(1,fx0,fy0)
!
!-----------------------------------------------------------------------
!
!  Establish coordinate for scalar forecast fields.
!
!-----------------------------------------------------------------------
!
    DO i=1,nx-1
      xs(i)=0.5*(x(i)+x(i+1))
    END DO
    xs(nx)=2.*xs(nx-1)-xs(nx-2)
    DO j=1,ny-1
      ys(j)=0.5*(y(j)+y(j+1))
    END DO
    ys(ny)=2.*ys(ny-1)-ys(ny-2)
    CALL xytoll(nx,ny,xs,ys,lat,lon)
    DO k=1,nz-1
      DO j=1,ny
        DO i=1,nx
          zps(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Get verification data file info including the name of the input data set.
!
!-----------------------------------------------------------------------
!
    READ(5,vrftn_fn,END=100)
    WRITE(6,'(a)')'Namelist block vrftn_fn sucessfully read in.' 
 
    lengbf=LEN_trim(vgrdbasfn)
    WRITE(6,'(/a,a)')' The grid/base name is ', vgrdbasfn(1:lengbf)

    lenfil=LEN_trim(vfilename)
    WRITE(6,'(/a,a)')' The data set name is ', vfilename(1:lenfil)
!
!-----------------------------------------------------------------------
!
!  Set the gridread parameter to 0 so that the verification
!  grid/base file will be read.
!
!-----------------------------------------------------------------------
!
    CALL setgbrd (0)
!
!-----------------------------------------------------------------------
!
!  Read in the verification data.
!
!-----------------------------------------------------------------------
!
    CALL dtaread(vnx,vny,vnz,nstyps,                                    &
                 vhinfmt,nch,vgrdbasfn(1:lengbf),lengbf,                &
                 vfilename(1:lenfil),lenfil,time,                       &
                 vx,vy,vz,vzp, vuprt ,vvprt ,vwprt ,vptprt, vpprt ,     &
                 vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,       &
                 vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar,   &
                 vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg,          &
                 vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,vsnowdpth,        &
                 vraing,vrainc,vprcrate,                                &
                 vradfrc,vradsw,vrnflx,                                 &
                 vusflx,vvsflx,vptsflx,vqvsflx,                         &
                 ireturn, vtem1,vtem2,vtem3)

    ivproj=mapproj
    vscale=sclfct
    vlatnot(1)=trulat1
    vlatnot(2)=trulat2
    vtrulon=trulon
    vdx=vx(3)-vx(2)
    vdy=vy(3)-vy(2)
    vctrlat=ctrlat
    vctrlon=ctrlon
    CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
    CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr)
    vx0=xctr-vdx*((vnx-3)/2)
    vy0=yctr-vdy*((vny-3)/2)
    CALL setorig(1,vx0,vy0)
!
!-----------------------------------------------------------------------
!
!  Establish coordinate for scalar verification fields.
!
!-----------------------------------------------------------------------
!
    DO i=1,vnx-1
      vxs(i)=0.5*(vx(i)+vx(i+1))
    END DO
    vxs(vnx)=2.*vxs(vnx-1)-vxs(vnx-2)
    DO j=1,vny-1
      vys(j)=0.5*(vy(j)+vy(j+1))
    END DO
    vys(vny)=2.*vys(vny-1)-vys(vny-2)
    DO k=1,vnz-1
      DO j=1,vny
        DO i=1,vnx
          vzps(i,j,k)=0.5*(vzp(i,j,k)+vzp(i,j,k+1))
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Find location of scalar forecast fields in verification grid.
!
!-----------------------------------------------------------------------
!
    IF(fx0 == vx0 .AND. fy0 == vy0 .AND.                                &
          flatnot(1) == vlatnot(1) .AND. flatnot(2) == vlatnot(2) .AND. &
          ftrulon == vtrulon .AND. ifproj == ivproj .AND.               &
          fscale == vscale ) THEN
      comcoord=.true.
      WRITE(6,'(//a//)') ' Grids share a common coordinate system'
      DO j=1,ny
        DO i=1,nx
          x2d(i,j)=xs(i)
        END DO
      END DO
      DO j=1,ny
        DO i=1,nx
          y2d(i,j)=ys(j)
        END DO
      END DO
    ELSE
      comcoord=.false.
      WRITE(6,'(//a,a//)') ' Grid coordinate systems differ',           &
                         ' Will convert coordinates via lat,lon.'
      CALL lltoxy(nx,ny,lat,lon,x2d,y2d)
    END IF
!
!-----------------------------------------------------------------------
!
!  Interpolate verification scalars to forecast grid.
!
!-----------------------------------------------------------------------
!

!EMK Special for same grids
     IF ((comcoord) .AND. (nx == vnx) .AND. (ny == vny) .AND. &
         (nz == vnz)) THEN
       comcoord2 = .TRUE.
       DO k = 1,nz
         IF (z(k) /= vz(k)) comcoord2 = .FALSE.
       END DO
       
       IF (comcoord2) THEN
         aptprt = vptprt
         apprt = vpprt
         aqvprt = vqvprt
         aqc = vqc
         aqr = vqr
         aqi = vqi
         aqs = vqs
         aqh = vqh
         atke = vtke
         akmh = vkmh
         akmv = vkmv
         aptbar = vptbar
         apbar = vpbar
         arhobar = vrhobar
         aqvbar = vqvbar
         atsfc = vtsfc
         atsoil = vtsoil
         awetdp = vwetdp
         awetcanp = vwetcanp
         awetsfc = vwetsfc
         araing = vraing
         arainc = vrainc
         aprcrate = vprcrate
         aradfrc = vradfrc
         aradsw = vradsw
         arnflx = vrnflx
         ausflx = vusflx
         avsflx = vvsflx
         aptsflx = vptsflx
         aqvsflx = vqvsflx
         awbar = vwbar
         awprt = vwprt
         avbar = vvbar
         avprt = vvprt
         aubar = vubar
         auprt = vuprt
       
         GO TO 1000 ! Skip interpolation
      END IF
    END IF
!EMK Special for same grids

    CALL setdxdy(vnx,vny,                                               &
                 1,vnx-1,1,vny-1,                                       &
                 vxs,vys,dxfld,dyfld,rdxfld,rdyfld)

    CALL intsclrs(nx, ny, nz, vnx, vny, vnz,                            &
                  1, nx-1,1, ny-1,1, nz-1,                              &
                  1,vnx-1,1,vny-1,1,vnz-1,                              &
                  iorder,                                               &
                  x2d, y2d, zps,vxs,vys,vzps,                           &
                  vptprt, vpprt,                                        &
                  vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,      &
                  vptbar, vpbar, vrhobar, vqvbar,                       &
                  vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,                 &
                  vraing,vrainc,vprcrate,                               &
                  vradfrc,vradsw,vrnflx,                                &
                  vusflx,vvsflx,vptsflx,vqvsflx,                        &
                  aptprt, apprt,                                        &
                  aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,      &
                  aptbar, apbar, arhobar, aqvbar,                       &
                  atsfc,atsoil,awetsfc,awetdp,awetcanp,                 &
                  araing,arainc,aprcrate,                               &
                  aradfrc,aradsw,arnflx,                                &
                  ausflx,avsflx,aptsflx,aqvsflx,                        &
                  iloc,jloc,zpver,dxfld,dyfld,rdxfld,rdyfld,            &
                  vtem1,vtem2,vtem3,                                    &
                  ireturn )
!
!-----------------------------------------------------------------------
!
!  Interpolate verification w to forecast grid.
!
!-----------------------------------------------------------------------
!
    CALL intonef(nx, ny, nz, vnx, vny, vnz,                             &
                 1,nx-1,1,ny-1,1,nz,                                    &
                 1,vnx-1,1,vny-1,1,vnz,                                 &
                 iorder,                                                &
                 x2d, y2d, zp, vxs,vys,vzp,                             &
                 vwprt, vwbar, awprt, awbar,                            &
                 iloc,jloc,zpver,dxfld,dyfld,rdxfld,rdyfld,             &
                 vtem1,vtem2,vtem3,                                     &
                 ireturn )
!
!-----------------------------------------------------------------------
!
!  Find location of u forecast field in verification grid.
!
!-----------------------------------------------------------------------
!
    IF(comcoord) THEN
      DO j=1,ny
        DO i=1,nx
          x2d(i,j)=x(i)
        END DO
      END DO
    ELSE
      CALL setmapr(ifproj,fscale,flatnot,ftrulon)
      CALL setorig(1,fx0,fy0)
      CALL xytoll(nx,ny,x,ys,lat,lon)
      CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
      CALL setorig(1,vx0,vy0)
      CALL lltoxy(nx,ny,lat,lon,x2d,y2d)
    END IF
!
!-----------------------------------------------------------------------
!
!  Interpolate verification u to forecast grid.
!
!-----------------------------------------------------------------------
!
    CALL setdxdy(vnx,vny,                                               &
                 1,vnx,1,vny-1,                                         &
                 vx,vys,dxfld,dyfld,rdxfld,rdyfld)
    CALL intonef(nx, ny, nz, vnx, vny, vnz,                             &
                 1,nx,1,ny-1,1,nz-1,                                    &
                 1,vnx,1,vny-1,1,vnz-1,                                 &
                 iorder,                                                &
                 x2d, y2d, zp, vxs,vys,vzp,                             &
                 vuprt, vubar, auprt, aubar,                            &
                 iloc,jloc,zpver,dxfld,dyfld,rdxfld,rdyfld,             &
                 vtem1,vtem2,vtem3,                                     &
                 ireturn )
!
!-----------------------------------------------------------------------
!
!  Find location of v forecast field in verification grid.
!
!-----------------------------------------------------------------------
!
    IF(comcoord) THEN
      DO j=1,ny
        DO i=1,nx
          x2d(i,j)=xs(i)
        END DO
      END DO
      DO j=1,ny
        DO i=1,nx
          y2d(i,j)=y(j)
        END DO
      END DO
    ELSE
      CALL setmapr(ifproj,fscale,flatnot,ftrulon)
      CALL setorig(1,fx0,fy0)
      CALL xytoll(nx,ny,xs,y,lat,lon)
      CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
      CALL setorig(1,vx0,vy0)
      CALL lltoxy(nx,ny,lat,lon,x2d,y2d)
    END IF
!
!-----------------------------------------------------------------------
!
!  Interpolate verification v to forecast grid.
!
!-----------------------------------------------------------------------
!
    CALL setdxdy(vnx,vny,                                               &
                 1,vnx-1,1,vny,                                         &
                 vxs,vy,dxfld,dyfld,rdxfld,rdyfld)
    CALL intonef(nx, ny, nz, vnx, vny, vnz,                             &
                 1,nx-1,1,ny,1,nz-1,                                    &
                 1,vnx-1,1,vny,1,vnz-1,                                 &
                 iorder,                                                &
                 x2d, y2d, zp, vxs,vys,vzp,                             &
                 vvprt, vvbar, avprt, avbar,                            &
                 iloc,jloc,zpver,dxfld,dyfld,rdxfld,rdyfld,             &
                 vtem1,vtem2,vtem3,                                     &
                 ireturn )
!
!-----------------------------------------------------------------------
!
!  Find   difference = forecast - verification
!
!  To reduce memory requirements, the difference fields are
!  written to the same arrays as the interpolated fields.
!
!-----------------------------------------------------------------------
!
    1000 CONTINUE

    CALL diffield(nx,ny,nz,                                             &
                  uprt, vprt, wprt, ptprt, pprt,                        &
                  qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,               &
                  ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,         &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,                      &
                  raing,rainc,prcrate,                                  &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx,                            &
                  auprt, avprt, awprt, aptprt, apprt,                   &
                  aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,      &
                  aubar, avbar, awbar, aptbar, apbar, arhobar, aqvbar,  &
                  atsfc,atsoil,awetsfc,awetdp,awetcanp,                 &
                  araing,arainc,aprcrate,                               &
                  aradfrc,aradsw,arnflx,                                &
                  ausflx,avsflx,aptsflx,aqvsflx,                        &
                  auprt, avprt, awprt, aptprt, apprt,                   &
                  aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,      &
                  tem1,atsfc,atsoil,awetsfc,awetdp,awetcanp,            &
                  araing,arainc,aprcrate,                               &
                  aradfrc,aradsw,arnflx,                                &
                  ausflx,avsflx,aptsflx,aqvsflx,                        &
                  ireturn )
!
!-----------------------------------------------------------------------
!
!  Set output variables to forecast coordinates
!
!-----------------------------------------------------------------------
!
    curtim=time
    mapproj=ifproj
    sclfct=fscale
    trulat1=flatnot(1)
    trulat2=flatnot(2)
    trulon=ftrulon
    ctrlat=fctrlat
    ctrlon=fctrlon
!
!-----------------------------------------------------------------------
!
!  Get output info
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/a/)')' Please tell the program about the OUTPUT file'
!
!-----------------------------------------------------------------------
!
!  Get runname to use for output data.
!
!-----------------------------------------------------------------------
!
    READ(5,output,END=100)
    WRITE(6,'(/5x,a,a)') 'The output run name is: ', runnmin

    runname=runnmin
!
!-----------------------------------------------------------------------
!
!  Find out the number of characters to be used to construct file
!  names.
!
!-----------------------------------------------------------------------
!
    CALL gtlfnkey( runname, lfnkey )
!
!-----------------------------------------------------------------------
!
!  Find out the number of characters to be used to construct file
!  names.
!
!-----------------------------------------------------------------------
!
    CALL gtlfnkey( runname, lfnkey )
!
    WRITE(6,'(a)')' Please input the data format flag value 1/2/3/4 '
!
!-----------------------------------------------------------------------
!
!  Set control parameters for
!  grid, base state, moisture, and ice variable dumping.
!
!-----------------------------------------------------------------------
!
    varout=1
    WRITE(6,'(a/a/a)')                                                  &
        ' Will it contain any grid information? (1 or 0)',              &
        ' If it will not, grid and base information will be dumped',    &
        ' to a separate file (filename input later as grdbasfn).'

    WRITE(6,'(a)')                                                      &
        ' Will it contain any base state data? (1 or 0)'

    WRITE(6,'(2(/5x,a)/)')                                              &
         'Do you want to write u, v, w, ptprt and pprt arrays?',        &
         '(select 0 or 1)'

    WRITE(6,'(a)')                                                      &
             ' Write moisture fields to the output file?(1 or 0)'

    WRITE(6,'(a)')' Write ice fields to the output file?(1 or 0)'

    WRITE(6,'(a)')                                                      &
             ' Write turbulence fields to the output file?(1 or 0)'

    WRITE(6,'(a)')                                                      &
        ' Write surface (soil) fields to the output file?(1 or 0)'

    WRITE(6,'(a)')' Write rain fields to the output file?(1 or 0)'

    WRITE(6,'(a)')' Write precipitation rate to output file?(1 or 0)'

    WRITE(6,'(a)')' Write radiation arrays to output file?(1 or 0)'

    WRITE(6,'(a)')' Write the surface fluxes to output file?(1 or 0)'

    WRITE(6,'(/5x,a/)')                                                 &
         'Do you want to compress the output data? (select 0 or 1)'


    CALL gtbasfn(runname(1:lfnkey),'./',2,hdmpfmt,mgrid,nestgrd,        &
                 grdbasfn, lengbf)

    WRITE(6,'(/1x,a,a)')                                                &
        'Output grid/base state file is ', grdbasfn(1:lengbf)

    nchdmp = 80
    grdbas = 1      ! Dump out grd and base state arrays only

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          auprt(i,j,k)=aubar(i,j,k)+auprt(i,j,k)
          avprt(i,j,k)=avbar(i,j,k)+avprt(i,j,k)
          awprt(i,j,k)=awbar(i,j,k)+awprt(i,j,k)
          aqvprt(i,j,k)=aqvbar(i,j,k)+aqvprt(i,j,k)
        END DO
      END DO
    END DO

    CALL dtadump(nx,ny,nz,nstyps,                                       &
                 hdmpfmt,nchdmp,grdbasfn(1:lengbf),grdbas,filcmprs,     &
                 auprt,avprt,awprt,aptprt,apprt,                        &
                 aqvprt,aqc,aqr,aqi,aqs,aqh,atke,akmh,akmv,             &
                 ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,                &
                 x,y,z,zp,hterain,j1,j2,j3,                             &
                 asoiltyp,astypfrct,avegtyp,alai,aroufns,aveg,          &
                 atsfc,atsoil,awetsfc,awetdp,awetcanp,asnowdpth,        &
                 araing,arainc,aprcrate,                                &
                 aradfrc,aradsw,arnflx,                                 &
                 ausflx,avsflx,aptsflx,aqvsflx,                         &
                 tem1,tem2,tem3)

!
!-----------------------------------------------------------------------
!
!  Find a unique name hdmpfn(1:ldmpf) for history dump data set
!  at time 'curtim'.
!
!-----------------------------------------------------------------------
!
    CALL gtdmpfn(runname(1:lfnkey),'./',2,                              &
                 curtim,hdmpfmt,                                        &
                 mgrid,nestgrd, hdmpfn, ldmpf)

    WRITE(6,'(/1x,a,f10.0,a,a)')                                        &
        'Output file at time ',curtim,' (s) is ', hdmpfn(1:ldmpf)

    grdbas = 0      ! Not just dump out grd and base state arrays only

    CALL dtadump(nx,ny,nz,nstyps,                                       &
                 hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,filcmprs,        &
                 auprt,avprt,awprt,aptprt,apprt,                        &
                 aqvprt,aqc,aqr,aqi,aqs,aqh,atke,akmh,akmv,             &
                 ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,                &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 asoiltyp,astypfrct,avegtyp,alai,aroufns,aveg,          &
                 atsfc,atsoil,awetsfc,awetdp,awetcanp,asnowdpth,        &
                 araing,arainc,aprcrate,                                &
                 aradfrc,aradsw,arnflx,                                 &
                 ausflx,avsflx,aptsflx,aqvsflx,                         &
                 tem1,tem2,tem3)

  END IF

  GOTO 101
  
  100 CONTINUE
  
  WRITE(6,'(a)')'Namelist block READ in error. Then program will terminated.'
  
  101 CONTINUE

  STOP

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


SUBROUTINE diffield(nx,ny,nz,                                           & 1,28
           uprt, vprt, wprt, ptprt, pprt,                               &
           qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                      &
           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,                &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,                             &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           auprt, avprt, awprt, aptprt, apprt,                          &
           aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,             &
           aubar, avbar, awbar, aptbar, apbar, arhobar, aqvbar,         &
           atsfc,atsoil,awetsfc,awetdp,awetcanp,                        &
           araing,arainc,aprcrate,                                      &
           aradfrc,aradsw,arnflx,                                       &
           ausflx,avsflx,aptsflx,aqvsflx,                               &
           duprt, dvprt, dwprt, dptprt, dpprt,                          &
           dqvprt, dqc, dqr, dqi, dqs, dqh, dtke,dkmh,dkmv,tem1,        &
           dtsfc,dtsoil,dwetsfc,dwetdp,dwetcanp,                        &
           draing,drainc,dprcrate,                                      &
           dradfrc,dradsw,drnflx,                                       &
           dusflx,dvsflx,dptsflx,dqvsflx,                               &
           ireturn)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Subtract the forecast fields from the interpolated verification
!  fields (names beginning with "a") and output to the difference
!  fields (names beginning with "d").  The input and difference
!  fields may share the same storage location.  For this subroutine
!  it is assumed the forecast and corresponding verification
!  data are at the same physical location, however, the physical
!  location may differ between variables. That is uprt and auprt
!  are at the same location, but that may differ from pprt and apprt.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster  Ou School of Meteorology. April 1992
!
!  MODIFICATION HISTORY:
!   14 May 1992  (KB) changed from arps2.5 to arps3.0
!   03 Aug 1992  (KB) updated to account for changes in arps3.0
!
!   09/07/1995  (KB)
!   Added differencing of surface (soil) fields.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!    nx,ny,nz Array dimensions for forecast field.
!
!    FORECAST FIELDS:
!
!    uprt     perturbation x component of velocity (m/s)
!    vprt     perturbation y component of velocity (m/s)
!    wprt     perturbation vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    kmh      Horizontal turbulent mixing coefficient (m**2/s)
!    kmv      Vertical turbulent mixing coefficient (m**2/s)
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    INTERPOLATED VERIFICATION FIELDS:
!
!    auprt     perturbation x component of velocity (m/s)
!    avprt     perturbation y component of velocity (m/s)
!    awprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    aptprt    perturbation potential temperature (K)
!    apprt     perturbation pressure (Pascal)
!
!    aqvprt    perturbation water vapor mixing ratio (kg/kg)
!    aqc       Cloud water mixing ratio (kg/kg)
!    aqr       Rainwater mixing ratio (kg/kg)
!    aqi       Cloud ice mixing ratio (kg/kg)
!    aqs       Snow mixing ratio (kg/kg)
!    aqh       Hail mixing ratio (kg/kg)
!
!    aubar     Base state x velocity component (m/s)
!    avbar     Base state y velocity component (m/s)
!    awbar     Base state z velocity component (m/s)
!    aptbar    Base state potential temperature (K)
!    apbar     Base state pressure (Pascal)
!    arhobar   Base state density (kg/m**3)
!    aqvbar    Base state water vapor mixing ratio (kg/kg)
!
!    atsfc     Temperature at surface (K)
!    atsoil    Deep soil temperature (K)
!    awetsfc   Surface soil moisture
!    awetdp    Deep soil moisture
!    awetcanp  Canopy water amount
!
!    araing    Grid supersaturation rain
!    arainc    Cumulus convective rain
!    aprcrate  Precipitation rates
!
!    aradfrc   Radiation forcing (K/s)
!    aradsw    Solar radiation reaching the surface
!    arnflx    Net radiation flux absorbed by surface
!
!    ausflx    Surface flux of u-momentum (kg/(m*s**2))
!    avsflx    Surface flux of v-momentum (kg/(m*s**2))
!    aptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    aqvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT :
!
!    DIFFERENCE FIELDS (may share storage with forecast fields
!                       or interpolated fields in calling program):
!
!    duprt     perturbation x component of velocity (m/s)
!    dvprt     perturbation y component of velocity (m/s)
!    dwprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    dptprt    perturbation potential temperature (K)
!    dpprt     perturbation pressure (Pascal)
!
!    dqvprt    perturbation water vapor mixing ratio (kg/kg)
!    dqc       Cloud water mixing ratio (kg/kg)
!    dqr       Rainwater mixing ratio (kg/kg)
!    dqi       Cloud ice mixing ratio (kg/kg)
!    dqs       Snow mixing ratio (kg/kg)
!    dqh       Hail mixing ratio (kg/kg)
!
!    dtsfc     Temperature at surface (K)
!    dtsoil    Deep soil temperature (K)
!    dwetsfc   Surface soil moisture
!    dwetdp    Deep soil moisture
!    dwetcanp  Canopy water amount
!
!    draing    Grid supersaturation rain
!    drainc    Cumulus convective rain
!    dprcrate  Precipitation rates
!
!    dradfrc   Radiation forcing (K/s)
!    dradsw    Solar radiation reaching the surface
!    drnflx    Net radiation flux absorbed by surface
!
!    dusflx    Surface flux of u-momentum (kg/(m*s**2))
!    dvsflx    Surface flux of v-momentum (kg/(m*s**2))
!    dptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    dqvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    tem1      Work array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! 3 dimensions of array
!
!-----------------------------------------------------------------------
!
!  Model Arrays
!
!-----------------------------------------------------------------------
!
  REAL :: uprt   (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: vprt   (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: wprt   (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: ptprt  (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: qvprt  (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: qc     (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: qr     (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: qi     (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs     (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: qh     (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

  REAL :: tke    (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: kmh    (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv    (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

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

  REAL :: tsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: tsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: wetsfc(nx,ny)       ! Surface soil moisture
  REAL :: wetdp (nx,ny)       ! Deep soil moisture
  REAL :: wetcanp(nx,ny)      ! Canopy water amount

  REAL :: raing (nx,ny)       ! Grid supersaturation rain
  REAL :: rainc (nx,ny)       ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Verification data interpolated to model grid
!
!-----------------------------------------------------------------------
!
  REAL :: auprt  (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: avprt  (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: awprt  (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: aptprt (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: apprt  (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: aqvprt (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: aqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: aqr    (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: aqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: aqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: aqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

  REAL :: atke   (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: akmh   (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: akmv   (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: aubar  (nx,ny,nz)    ! Base state u-velocity (m/s)
  REAL :: avbar  (nx,ny,nz)    ! Base state v-velocity (m/s)
  REAL :: awbar  (nx,ny,nz)    ! Base state w-velocity (m/s)
  REAL :: aptbar (nx,ny,nz)    ! Base state potential temperature (K)
  REAL :: arhobar(nx,ny,nz)    ! Base state density (kg/m**3)
  REAL :: apbar  (nx,ny,nz)    ! Base state pressure (Pascal)
  REAL :: aqvbar (nx,ny,nz)    ! Base state water vapor specific humidity

  REAL :: atsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: atsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: awetsfc(nx,ny)       ! Surface soil moisture
  REAL :: awetdp (nx,ny)       ! Deep soil moisture
  REAL :: awetcanp(nx,ny)      ! Canopy water amount

  REAL :: araing (nx,ny)       ! Grid supersaturation rain
  REAL :: arainc (nx,ny)       ! Cumulus convective rain
  REAL :: aprcrate(nx,ny,4)    ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: aradfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: aradsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: arnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: ausflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: avsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: aptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: aqvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Difference arrays
!
!-----------------------------------------------------------------------
!
  REAL :: duprt  (nx,ny,nz)    ! perturbation x component of velocity (m/s)
  REAL :: dvprt  (nx,ny,nz)    ! perturbation y component of velocity (m/s)
  REAL :: dwprt  (nx,ny,nz)    ! perturbation vertical component of
                               ! velocity in Cartesian coordinates (m/s)
  REAL :: dptprt (nx,ny,nz)    ! perturbation potential temperature (K)
  REAL :: dpprt  (nx,ny,nz)    ! perturbation pressure (Pascal)
  REAL :: dqvprt (nx,ny,nz)    ! perturbation water vapor mixing ratio (kg/kg)
  REAL :: dqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: dqr    (nx,ny,nz)    ! Rainwater mixing ratio (kg/kg)
  REAL :: dqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: dqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: dqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)
  REAL :: dtke   (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: dkmh   (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: dkmv   (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: dtsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: dtsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: dwetsfc(nx,ny)       ! Surface soil moisture
  REAL :: dwetdp (nx,ny)       ! Deep soil moisture
  REAL :: dwetcanp(nx,ny)      ! Canopy water amount

  REAL :: draing (nx,ny)       ! Grid supersaturation rain
  REAL :: drainc (nx,ny)       ! Cumulus convective rain
  REAL :: dprcrate(nx,ny,4)    ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: dradfrc(nx,ny,nz)    ! Radiation forcing (K/s)
  REAL :: dradsw (nx,ny)       ! Solar radiation reaching the surface
  REAL :: drnflx (nx,ny)       ! Net radiation flux absorbed by surface

  REAL :: dusflx (nx,ny)       ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: dvsflx (nx,ny)       ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: dptsflx(nx,ny)       ! Surface heat flux (K*kg/(m*s**2))
  REAL :: dqvsflx(nx,ny)       ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1   (nx,ny,nz)    ! A work array
  INTEGER :: ireturn, i,j,k
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: is,js,ks,ie,je,ke
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  is=1
  js=1
  ks=1
  ie=nx-1
  je=ny-1
  ke=nz-1

!-----------------------------------------------------------------------
!
!  Scalars
!
!-----------------------------------------------------------------------

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem1(i,j,k)=0.0
      END DO
    END DO
  END DO

  PRINT *, ' ptprt: '
  CALL subtr(nx,ny,nz, ptprt,ptbar,aptprt,aptbar,dptprt,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' pprt: '
  CALL subtr(nx,ny,nz,  pprt, pbar, apprt, apbar, dpprt,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qvprt: '
  CALL subtr(nx,ny,nz, qvprt,qvbar,aqvprt,aqvbar,dqvprt,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qc: '
  CALL subtr(nx,ny,nz,    qc, tem1,   aqc,  tem1,   dqc,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qr: '
  CALL subtr(nx,ny,nz,    qr, tem1,   aqr,  tem1,   dqr,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qi: '
  CALL subtr(nx,ny,nz,    qi, tem1,   aqi,  tem1,   dqi,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qs: '
  CALL subtr(nx,ny,nz,    qs, tem1,   aqs,  tem1,   dqs,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' qh: '
  CALL subtr(nx,ny,nz,    qh, tem1,   aqh,  tem1,   dqh,                &
             is,js,ks,ie,je,ke)
  PRINT *, ' tke: '
  CALL subtr(nx,ny,nz,    tke, tem1,   atke,  tem1,   dtke,             &
             is,js,ks,ie,je,ke)
  PRINT *, ' kmh: '
  CALL subtr(nx,ny,nz,    kmh, tem1,   akmh,  tem1,   dkmh,             &
             is,js,ks,ie,je,ke)
  PRINT *, ' kmv: '
  CALL subtr(nx,ny,nz,    kmv, tem1,   akmv,  tem1,   dkmv,             &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  u wind components
!
!-----------------------------------------------------------------------

  ie=nx
  PRINT *, ' uprt: '
  CALL subtr(nx,ny,nz,uprt,ubar,auprt,aubar,duprt,                      &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  v wind components
!
!-----------------------------------------------------------------------

  ie=nx-1
  je=ny
  PRINT *, ' vprt: '
  CALL subtr(nx,ny,nz,vprt,vbar,avprt,avbar,dvprt,                      &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  w wind components
!
!-----------------------------------------------------------------------

  je=ny-1
  ke=nz
  PRINT *, ' wprt: '
  CALL subtr(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,                       &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  2-d surface (soil) variables
!
!-----------------------------------------------------------------------

  ie=nx-1
  je=ny-1
  ks=1
  ke=1
!
  PRINT *, ' tsfc:'
  CALL subtr(nx,ny,1,   tsfc,tem1,   atsfc,tem1,   dtsfc,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' tsoil:'
  CALL subtr(nx,ny,1,  tsoil,tem1,  atsoil,tem1,  dtsoil,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetsfc:'
  CALL subtr(nx,ny,1, wetsfc,tem1, awetsfc,tem1, dwetsfc,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetdp:'
  CALL subtr(nx,ny,1,  wetdp,tem1,  awetdp,tem1,  dwetdp,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetcanp:'
  CALL subtr(nx,ny,1,wetcanp,tem1,awetcanp,tem1,dwetcanp,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' raing:'
  CALL subtr(nx,ny,1,  raing,tem1,  araing,tem1,  draing,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' rainc:'
  CALL subtr(nx,ny,1,  rainc,tem1,  arainc,tem1,  drainc,               &
             is,js,ks,ie,je,ke)
  PRINT *, ' prcrate1:'
  CALL subtr(nx,ny,1,  prcrate(1,1,1),tem1, aprcrate(1,1,1),tem1,       &
             dprcrate(1,1,1), is,js,ks,ie,je,ke)
  PRINT *, ' prcrate2:'
  CALL subtr(nx,ny,1,  prcrate(1,1,2),tem1, aprcrate(1,1,2),tem1,       &
             dprcrate(1,1,2), is,js,ks,ie,je,ke)
  PRINT *, ' prcrate3:'
  CALL subtr(nx,ny,1,  prcrate(1,1,3),tem1, aprcrate(1,1,3),tem1,       &
             dprcrate(1,1,3), is,js,ks,ie,je,ke)
  PRINT *, ' prcrate4:'
  CALL subtr(nx,ny,1,  prcrate(1,1,4),tem1, aprcrate(1,1,4),tem1,       &
             dprcrate(1,1,4), is,js,ks,ie,je,ke)

  PRINT *, ' radfrc:'
  CALL subtr(nx,ny,nz, radfrc,tem1, aradfrc,tem1,                       &
             dradfrc, is,js,ks,ie,je,ke)
  PRINT *, ' radsw:'
  CALL subtr(nx,ny,1, radsw,tem1, aradsw,tem1,                          &
             dradsw, is,js,ks,ie,je,ke)
  PRINT *, ' rnflx:'
  CALL subtr(nx,ny,1, rnflx,tem1, arnflx,tem1,                          &
             drnflx, is,js,ks,ie,je,ke)

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


SUBROUTINE subtr(nx,ny,nz, a,abar,b,bbar,c,                             & 28
           istr,jstr,kstr,iend,jend,kend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Subtracts 2 three-dimensional arrays, represented by
!  means plus perturbations.
!
!  AUTHOR: Keith Brewster  OU School of Meteorology.  Feb 1992
!
!  MODIFICATION HISTORY:
!   11 Aug 1992  (KB) changed from arps2.5 to arps3.0
!
!   7 May 2002 (Eric Kemp)
!   Minor change to list all non-zero differences.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    a        perturbation data array
!    abar     mean data array
!    b        perturbation data array to subtract from a
!    bbar     mean data array to subtract from a
!
!  OUTPUT:
!    c        difference array a-b
!             (may share storage in calling program with array a or b)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! 3 dimensions of array

  REAL :: a   (nx,ny,nz)       ! data array
  REAL :: abar(nx,ny,nz)       ! base state of data array a
  REAL :: b   (nx,ny,nz)       ! data array to subtract from a
  REAL :: bbar(nx,ny,nz)       ! base state of data arrya b
  REAL :: c   (nx,ny,nz)       ! difference array a-b
  INTEGER :: istr,jstr,kstr
  INTEGER :: iend,jend,kend
  INTEGER :: i,j,k,imid,jmid,kmid
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=nint(0.5*FLOAT(istr+iend))
  jmid=nint(0.5*FLOAT(jstr+jend))
  kmid=nint(0.5*FLOAT(kstr+kend))
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
  PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)),     &
                   ' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
!
!-----------------------------------------------------------------------
!
!  Subtraction
!
!-----------------------------------------------------------------------
!
  DO k=kstr,kend
    DO j=jstr,jend
      DO i=istr,iend
        c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k))

          IF (c(i,j,k) /= 0) THEN
            WRITE(6,*)'Non-zero difference in field! c = ',c(i,j,k), &
                      ' at i,j,k: ',i,j,k
          END IF 

      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
  PRINT *, '         c= ',c(imid,jmid,kmid)
  RETURN
END SUBROUTINE subtr
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE INTSCLRS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE intsclrs(nx,ny,nz,vnx,vny,vnz,                               & 1,34
           ibeg,iend,jbeg,jend,kbeg,kend,                               &
           ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                         &
           iorder,                                                      &
           xs2d,ys2d,zps, vxs,vys,vzps,                                 &
           vptprt, vpprt,                                               &
           vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,             &
           vptbar, vpbar, vrhobar, vqvbar,                              &
           vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,                        &
           vraing,vrainc,vprcrate,                                      &
           vradfrc,vradsw,vrnflx,                                       &
           vusflx,vvsflx,vptsflx,vqvsflx,                               &
           aptprt, apprt,                                               &
           aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,             &
           aptbar, apbar, arhobar, aqvbar,                              &
           atsfc,atsoil,awetsfc,awetdp,awetcanp,                        &
           araing,arainc,aprcrate,                                      &
           aradfrc,aradsw,arnflx,                                       &
           ausflx,avsflx,aptsflx,aqvsflx,                               &
           iloc,jloc,zpver,dxfld,dyfld,rdxfld,rdyfld,                   &
           slopey,alphay,betay,                                         &
           ireturn )
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!    Intfield interpolates scalars from a set of fields (the verification
!    fields, "verif") having Cartesian coordinates described by vx,vy,vzp
!    to a second set of fields described by cartesion coordinates x,y,zp.
!    It is assumed that x,y,zp and vx,vy,vzp are monotonically increasing
!    with increasing index.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster  OU School of Meteorology.  Feb 1992
!
!  MODIFICATION HISTORY:
!    12 Aug 1992  (KB) changed from arps2.5 to arps3.0
!    19 May 1993  (KB) changed from arps3.1 to arps3.2
!    24 May 1993  (KB) changed to special version for scalars only.
!
!     9 Sep 1995  (KB) added processing of sfc (soil) fields
!    26 Apr 1996  (KB) Version 2.0 -- Uses Gauss Forward routines for
!                      interpolation rather than piecewise linear.
!    07 Nov 1996  (KB) Replaced interpolation scheme.
!                      Reordered sequence of variables in call.
!    03 Nov 1999  (KB via Eric Kemp) Corrected dimensions of precip
!                      and flux v* (verif) variables that had been
!                      recently added.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    nx       Number of grid pts in the x-direction (east/west)
!    ny       Number of grid pts in the y-direction (north/south)
!    nz       Number of grid pts in the vertical
!
!    vnx      Number of verif grid points in the x-direction (east/west)
!    vny      Number of verif grid points in the y-direction (north/south)
!    vnz      Number of verif grid points in the vertical
!
!    ibeg,iend   Range of x index to do interpolation
!    jbeg,jend   Range of y index to do interpolation
!    kbeg,kend   Range of z index to do interpolation
!
!    ivbeg,ivend   Range of x index to use in verification array
!    jvbeg,jvend   Range of y index to use in verification array
!    kvbeg,kvend   Range of z index to use in verification array
!
!    iorder   Interpolation parameter.
!             iorder specifies the order of interpolation
!             1 = bi-linear
!             2 = bi-quadratic
!
!    xs2d     x coordinate of scalar grid points in physical space (m)
!    ys2d     y coordinate of scalar grid points in physical space (m)
!    zps      z coordinate of scalar grid points in physical space (m)
!
!    vxs      x coordinate of verif scalar grid points in physical space (m)
!    vys      y coordinate of verif scalar grid points in physical space (m)
!    vzps     z coordinate of verif scalar grid points in physical space (m)
!
!    vpt      Potential temperature
!    vpprt    Perturbation pressure  (Pascal)
!    vqv      Water vapor specific humidity  (kg/kg)
!    vqc      Cloud water mixing ratio  (kg/kg)
!    vqr      Rainwater mixing ratio  (kg/kg)
!    vqi      Cloud ice mixing ratio  (kg/kg)
!    vqs      Snow mixing ratio  (kg/kg)
!    vqh      Hail mixing ratio  (kg/kg)
!
!    vptbar    Base state potential temperature (K)
!    vpbar     Base state pressure (Pascal)
!    vrhobar   Base state density (kg/m**3)
!    vqvbar    Base state water vapor mixing ratio (kg/kg)
!
!    vtsfc     Temperature at surface (K)
!    vtsoil    Deep soil temperature (K)
!    vwetsfc   Surface soil moisture
!    vwetdp    Deep soil moisture
!    vwetcanp  Canopy water amount
!
!    vraing    Grid supersaturation rain
!    vrainc    Cumulus convective rain
!    vprcrate  Precipitation rates
!
!    vradfrc   Radiation forcing (K/s)
!    vradsw    Solar radiation reaching the surface
!    vrnflx    Net radiation flux absorbed by surface
!
!    vusflx    Surface flux of u-momentum (kg/(m*s**2))
!    vvsflx    Surface flux of v-momentum (kg/(m*s**2))
!    vptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    vqvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT:
!    apt      Interpolated potential temperature
!    apprt    Interpolated perturbation pressure  (Pascal)
!    aqv      Interpolated water vapor specific humidity  (kg/kg)
!    aqc      Interpolated cloud water mixing ratio  (kg/kg)
!    aqr      Interpolated rainwater mixing ratio  (kg/kg)
!    aqi      Interpolated cloud ice mixing ratio  (kg/kg)
!    aqs      Interpolated snow mixing ratio  (kg/kg)
!    aqh      Interpolated hail mixing ratio  (kg/kg)
!
!    aptbar   Interpolated base state potential temperature (K)
!    apbar    Interpolated base state pressure (Pascal)
!    arhobar  Interpolated base state density (kg/m**3)
!    aqvbar   Interpolated base state water vapor mixing ratio (kg/kg)
!
!    atsfc     Interpolated temperature at surface (K)
!    atsoil    Interpolated deep soil temperature (K)
!    awetsfc   Interpolated surface soil moisture
!    awetdp    Interpolated deep soil moisture
!    awetcanp  Interpolated canopy water amount
!
!    araing    Interpolated grid supersaturation rain
!    arainc    Interpolated cumulus convective rain
!    aprcrate  Precipitation rates
!
!    aradfrc   Radiation forcing (K/s)
!    aradsw    Solar radiation reaching the surface
!    arnflx    Net radiation flux absorbed by surface
!
!    ausflx    Surface flux of u-momentum (kg/(m*s**2))
!    avsflx    Surface flux of v-momentum (kg/(m*s**2))
!    aptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    aqvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    ireturn
!
!  WORK ARRAYS:
!
!    iloc     i-index of interpolation points in field to be interpolated
!    jloc     j-index of interpolation points in field to be interpolated
!    dxfld    Vector of delta-x (m) of field to be interpolated
!    dyfld    Vector of delta-y (m) of field to be interpolated
!    rdxfld   Vector of 1./delta-x (1/m) of field to be interpolated
!    rdyfld   Vector of 1./delta-y (1/m) of field to be interpolated
!
!    slopey   Piecewise linear df/dy
!    alphay   Coefficient of y-squared term in y quadratic interpolator
!    betay    Coefficient of y term in y quadratic interpolator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: vnx,vny,vnz
  INTEGER :: ibeg,iend,jbeg,jend,kbeg,kend
  INTEGER :: ivbeg,ivend,jvbeg,jvend,kvbeg,kvend
!
  REAL :: xs2d(nx,ny)
  REAL :: ys2d(nx,ny)
  REAL :: zps(nx,ny,nz)
!
  REAL :: vxs(vnx)
  REAL :: vys(vny)
  REAL :: vzps(vnx,vny,vnz)
!
  INTEGER :: iorder
!
!-----------------------------------------------------------------------
!
!  Original arrays (verification field)
!
!-----------------------------------------------------------------------
!
  REAL :: vptprt(vnx,vny,vnz)
  REAL :: vpprt (vnx,vny,vnz)
  REAL :: vqvprt(vnx,vny,vnz)
  REAL :: vqc   (vnx,vny,vnz)
  REAL :: vqr   (vnx,vny,vnz)
  REAL :: vqi   (vnx,vny,vnz)
  REAL :: vqs   (vnx,vny,vnz)
  REAL :: vqh   (vnx,vny,vnz)
  REAL :: vtke  (vnx,vny,vnz)
  REAL :: vkmh  (vnx,vny,vnz)
  REAL :: vkmv  (vnx,vny,vnz)

  REAL :: vptbar (vnx,vny,vnz)
  REAL :: vrhobar(vnx,vny,vnz)
  REAL :: vpbar  (vnx,vny,vnz)
  REAL :: vqvbar (vnx,vny,vnz)

  REAL :: vtsfc  (vnx,vny)     ! Temperature at surface (K)
  REAL :: vtsoil (vnx,vny)     ! Deep soil temperature (K)
  REAL :: vwetsfc(vnx,vny)     ! Surface soil moisture
  REAL :: vwetdp (vnx,vny)     ! Deep soil moisture
  REAL :: vwetcanp(vnx,vny)    ! Canopy water amount

  REAL :: vraing (vnx,vny)     ! Grid supersaturation rain
  REAL :: vrainc (vnx,vny)     ! Cumulus convective rain
  REAL :: vprcrate(vnx,vny,4)    ! precipitation rates (kg/(m**2*s))
                                 ! prcrate(1,1,1) = total precip. rate
                                 ! prcrate(1,1,2) = grid scale precip. rate
                                 ! prcrate(1,1,3) = cumulus precip. rate
                                 ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: vradfrc(vnx,vny,vnz)    ! Radiation forcing (K/s)
  REAL :: vradsw (vnx,vny)       ! Solar radiation reaching the surface
  REAL :: vrnflx (vnx,vny)       ! Net radiation flux absorbed by surface

  REAL :: vusflx (vnx,vny)       ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vvsflx (vnx,vny)       ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: vptsflx(vnx,vny)       ! Surface heat flux (K*kg/(m*s**2))
  REAL :: vqvsflx(vnx,vny)       ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Arrays interpolated to model grid
!
!-----------------------------------------------------------------------
!
  REAL :: aptprt(nx,ny,nz)
  REAL :: apprt (nx,ny,nz)
  REAL :: aqvprt(nx,ny,nz)
  REAL :: aqc   (nx,ny,nz)
  REAL :: aqr   (nx,ny,nz)
  REAL :: aqi   (nx,ny,nz)
  REAL :: aqs   (nx,ny,nz)
  REAL :: aqh   (nx,ny,nz)
  REAL :: atke  (nx,ny,nz)
  REAL :: akmh  (nx,ny,nz)
  REAL :: akmv  (nx,ny,nz)

  REAL :: aptbar (nx,ny,nz)
  REAL :: arhobar(nx,ny,nz)
  REAL :: apbar  (nx,ny,nz)
  REAL :: aqvbar (nx,ny,nz)

  REAL :: atsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: atsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: awetsfc(nx,ny)       ! Surface soil moisture
  REAL :: awetdp (nx,ny)       ! Deep soil moisture
  REAL :: awetcanp(nx,ny)      ! Canopy water amount

  REAL :: araing (nx,ny)       ! Grid supersaturation rain
  REAL :: arainc (nx,ny)       ! Cumulus convective rain
  REAL :: aprcrate(nx,ny,4)    ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: aradfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: aradsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: arnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: ausflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: avsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: aptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: aqvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Work arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: iloc(nx,ny)
  INTEGER :: jloc(nx,ny)
  REAL :: zpver(nx,ny,vnz)
!
  REAL :: dxfld(vnx)
  REAL :: dyfld(vny)
  REAL :: rdxfld(vnx)
  REAL :: rdyfld(vny)
  REAL :: slopey(vnx,vny,vnz)
  REAL :: alphay(vnx,vny,vnz)
  REAL :: betay(vnx,vny,vnz)
!
  INTEGER :: ireturn
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: k,korder
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Find i,j indices in verfication grid of each forecast point
!
!-----------------------------------------------------------------------
!
  CALL setijloc(nx,ny,vnx,vny,xs2d,ys2d,vxs,vys,iloc,jloc)
  CALL setdxdy(vnx,vny,                                                 &
               1,vnx-1,1,vny-1,                                         &
               vxs,vys,dxfld,dyfld,rdxfld,rdyfld)
!
!-----------------------------------------------------------------------
!
!  Interpolate 2-d fields
!
!-----------------------------------------------------------------------
!
  korder=MIN(iorder,2)
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vtsfc,vxs,vys,iloc,jloc,               &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                atsfc)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vtsoil,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                atsoil)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vwetdp,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                awetdp)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vwetsfc,vxs,vys,iloc,jloc,             &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                awetsfc)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vwetcanp,vxs,vys,iloc,jloc,            &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                awetcanp)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vraing,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                araing)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vrainc,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                arainc)
!
  DO k=1,4
    CALL fldint2d(nx,ny,vnx,vny,                                        &
                  ibeg,iend,jbeg,jend,                                  &
                  ivbeg,ivend,jvbeg,jvend,                              &
                  korder,xs2d,ys2d,vprcrate(1,1,k),                     &
                  vxs,vys,iloc,jloc,                                    &
                  dxfld,dyfld,rdxfld,rdyfld,                            &
                  slopey,alphay,betay,                                  &
                  aprcrate(1,1,k))
  END DO

  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vradsw,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aradsw)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vrnflx,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                arnflx)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vusflx,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                ausflx)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vvsflx,vxs,vys,iloc,jloc,              &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                avsflx)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vptsflx,vxs,vys,iloc,jloc,             &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aptsflx)
!
  CALL fldint2d(nx,ny,vnx,vny,                                          &
                ibeg,iend,jbeg,jend,                                    &
                ivbeg,ivend,jvbeg,jvend,                                &
                korder,xs2d,ys2d,vqvsflx,vxs,vys,iloc,jloc,             &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqvsflx)
!
!-----------------------------------------------------------------------
!
!  Create array of verification heights at
!  forecast x,y locations
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    CALL fldint2d(nx,ny,vnx,vny,                                        &
                  ibeg,iend,jbeg,jend,                                  &
                  ivbeg,ivend,jvbeg,jvend,                              &
                  korder,xs2d,ys2d,vzps(1,1,k),vxs,vys,iloc,jloc,       &
                  dxfld,dyfld,rdxfld,rdyfld,                            &
                  slopey,alphay,betay,                                  &
                  zpver(1,1,k))
  END DO
!
!-----------------------------------------------------------------------
!
!  Interpolate 3d scalar fields
!
!-----------------------------------------------------------------------
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vptprt,                            &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aptprt)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vpprt,                             &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                apprt)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqvprt,                            &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqvprt)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqc,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqc)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqr,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqr)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqi,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqi)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqs,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqs)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqh,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqh)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqi,                               &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqi)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vtke,                              &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                atke)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vkmh,                              &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                akmh)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vkmv,                              &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                akmv)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vptbar,                            &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aptbar)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vpbar,                             &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                apbar)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vrhobar,                           &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                arhobar)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vqvbar,                            &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aqvbar)
!
  CALL fldint3d(nx,ny,nz,vnx,vny,vnz,                                   &
                ibeg,iend,jbeg,jend,kbeg,kend,                          &
                ivbeg,ivend,jvbeg,jvend,kvbeg,kvend,                    &
                korder,xs2d,ys2d,zps,vradfrc,                           &
                vxs,vys,zpver,iloc,jloc,                                &
                dxfld,dyfld,rdxfld,rdyfld,                              &
                slopey,alphay,betay,                                    &
                aradfrc)
!
  RETURN
END SUBROUTINE intsclrs