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


PROGRAM arpsintrp,335
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This program interpolates gridded data from one ARPS grid to another.
!  It can be used to prepare data for running ARPS in a one-way nested
!  mode. It's exepcted to replace ARPSR2H in this capacity.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!
!  3/27/1997. Written based on ARPSR2H and ARPSCVT.
!
!  MODIFICATION HISTORY:
!
!  7/31/1997 (M. Xue)
!  Added options for specifying the terrain for the new grid, in
!  addition to interpolated terrain.
!
!  4/24/1998 (D. Weber)
!  Added hydrostatic extrapolation calculations for pprt and pbar
!  when the fine grid terrain is lower in elevation than the
!  coarse grid terrain. (DO LOOP 1423)  Note: this update is
!  functional for bglopt=2 ONLY.
!
!  4/24/1998 (D. Weber)
!  Added option for quadratic interpolation in the vertical.
!      intrpvopt = 1 for linear
!      intrpvopt = 2 for quadratic
!
!  4/1/1999 (M. Xue)
!
!  Rewrote major portions of the program. Special treatment given
!  to grid points below input grid ground. When output grid
!  terrain falls below ground, base-state variables are reconstructed
!  from averages on horizontal planes to ensure x or y independency
!  below ground. This is done only for bglopt=4, however.
!  The new code also runs much faster.
!
!  2000/04/05 (Gene Bassett)
!  Added bglopt=5, which is similar to bglopt=4 but uses the
!  method in ext2arps for constructing the base sounding.
!
!  2000/04/17 (Ming Xue)
!  Added an option that allows one to specify input history data
!  at a constant time interval.
!
!  2000/05/20 (Gene Bassett)
!  Converted to F90, creating allocation and main subroutines.
!
!  2000/07/28 (Ming Xue)
!  Converted to F90 free format. Use ALLOCATABLE instead of
!  POINTER allocation to avoid double memory usage.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!  2000/10/31 (Gene Bassett)
!  Refined handling of interpolating soil variables.
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!  2001/06/28 (Gene Bassett)
!  Added ntagopt, option to set surface winds in new grid to the surface
!  winds in the original grid for areas where new terrain is above the 
!  original terrain.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  REAL, allocatable :: u     (:,:,:)  ! Total u-velocity (m/s).
  REAL, allocatable :: v     (:,:,:)  ! Total v-velocity (m/s).
  REAL, allocatable :: w     (:,:,:)  ! Total w-velocity (m/s).
  REAL, allocatable :: ptprt (:,:,:)  ! Perturbation potential temperature
                                      ! from that of base state atmosphere (Kelvin).
  REAL, allocatable :: pprt  (:,:,:)  ! Perturbation pressure from that
                                      ! of base state atmosphere (Pascal).
  REAL, allocatable :: qv    (:,:,:)  ! Water vapor specific humidity (kg/kg).
  REAL, allocatable :: qc    (:,:,:)  ! Cloud water mixing ratio (kg/kg).
  REAL, allocatable :: qr    (:,:,:)  ! Rain water mixing ratio (kg/kg).
  REAL, allocatable :: qi    (:,:,:)  ! Cloud ice mixing ratio (kg/kg).
  REAL, allocatable :: qs    (:,:,:)  ! Snow mixing ratio (kg/kg).
  REAL, allocatable :: qh    (:,:,:)  ! Hail mixing ratio (kg/kg).
  REAL, allocatable :: 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 )

  INTEGER, allocatable :: soiltyp(:,:,:)   ! Soil type
  REAL, allocatable ::    stypfrct(:,:,:)  ! Fraction of 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 ::    ndvi   (:,:)     ! NDVI

  REAL, allocatable :: tsfc    (:,:,:)     ! Ground sfc. temperature (K)
  REAL, allocatable :: tsoil   (:,:,:)     ! Deep soil temperature (K)
  REAL, allocatable :: wetsfc  (:,:,:)     ! Surface soil moisture
  REAL, allocatable :: wetdp   (:,:,:)     ! Deep soil moisture
  REAL, allocatable :: wetcanp (:,:,:)     ! Canopy water amount
  REAL, allocatable :: snowdpth(:,:)       ! Snow depth (m)

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

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

  REAL, allocatable :: usflx (:,:)         ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, allocatable :: vsflx (:,:)         ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, allocatable :: ptsflx(:,:)         ! Surface heat flux (K*kg/(m*s**2))
  REAL, allocatable :: qvsflx(:,:)         ! Surface moisture flux (kg/(m**2*s))

  REAL, allocatable :: ubar  (:,:,:)       ! Base state u-velocity (m/s).
  REAL, allocatable :: vbar  (:,:,:)       ! Base state v-velocity (m/s).
  REAL, allocatable :: wbar  (:,:,:)       ! Base state w-velocity (m/s).
  REAL, allocatable :: ptbar (:,:,:)       ! Base state potential temperature (K)
  REAL, allocatable :: pbar  (:,:,:)       ! Base state pressure (Pascal).
  REAL, allocatable :: rhobar(:,:,:)       ! Base state air density (kg/m**3)
  REAL, allocatable :: qvbar (:,:,:)       ! Base state water vapor specific humidity
                                           ! (kg/kg).
  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 on the staggered grid.
  REAL, allocatable :: hterain(:,:)        ! The height of terrain.

  REAL, allocatable ::    uprt   (:,:,:)   ! Perturbation u-velocity (m/s)
  REAL, allocatable ::    vprt   (:,:,:)   ! Perturbation v-velocity (m/s)
  REAL, allocatable ::    qvprt  (:,:,:)   ! Perturbation water vapor specific
                                           ! humidity (kg/kg)

  REAL, allocatable :: tem1  (:,:,:)       ! Temporary array
  REAL, allocatable :: tem2  (:,:,:)       ! Temporary array
  REAL, allocatable :: tem3  (:,:,:)       ! Temporary array
  !wdt update
  REAL, allocatable :: tem4  (:,:,:)       ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Arrays on the new grid:
!
!-----------------------------------------------------------------------
!
  REAL, allocatable :: u1    (:,:,:)  ! Total u-velocity (m/s).
  REAL, allocatable :: v1    (:,:,:)  ! Total v-velocity (m/s).
  REAL, allocatable :: w1    (:,:,:)  ! Total w-velocity (m/s).
  REAL, allocatable :: ptprt1(:,:,:)  ! Perturbation potential temperature
                                      ! from that of base state atmosphere (Kelvin).
  REAL, allocatable :: pprt1 (:,:,:)  ! Perturbation pressure from that
                                      ! of base state atmosphere (Pascal).
  REAL, allocatable :: qv1   (:,:,:)  ! Water vapor specific humidity (kg/kg).
  REAL, allocatable :: qc1   (:,:,:)  ! Cloud water mixing ratio (kg/kg).
  REAL, allocatable :: qr1   (:,:,:)  ! Rain water mixing ratio (kg/kg).
  REAL, allocatable :: qi1   (:,:,:)  ! Cloud ice mixing ratio (kg/kg).
  REAL, allocatable :: qs1   (:,:,:)  ! Snow mixing ratio (kg/kg).
  REAL, allocatable :: qh1   (:,:,:)  ! Hail mixing ratio (kg/kg).
  REAL, allocatable :: tke1  (:,:,:)  ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, allocatable :: kmh1  (:,:,:)  ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, allocatable :: kmv1  (:,:,:)  ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )

  INTEGER, allocatable :: soiltyp1(:,:,:)   ! Soil type
  REAL, allocatable ::   stypfrct1(:,:,:)   ! Frction of soil type
  INTEGER, allocatable :: vegtyp1 (:,:)     ! Vegetation type
  REAL, allocatable ::    lai1    (:,:)     ! Leaf Area Index
  REAL, allocatable ::    roufns1 (:,:)     ! Surface roughness
  REAL, allocatable ::    veg1    (:,:)     ! Vegetation fraction

  REAL, allocatable :: tsfc1   (:,:,:)      ! Ground sfc. temperature (K)
  REAL, allocatable :: tsoil1  (:,:,:)      ! Deep soil temperature (K)
  REAL, allocatable :: wetsfc1 (:,:,:)      ! Surface soil moisture
  REAL, allocatable :: wetdp1  (:,:,:)      ! Deep soil moisture
  REAL, allocatable :: wetcanp1(:,:,:)      ! Canopy water amount
  REAL, allocatable :: snowdpth1(:,:)       ! Snow depth (m)

  REAL, allocatable :: raing1(:,:)          ! Grid supersaturation rain
  REAL, allocatable :: rainc1(:,:)          ! Cumulus convective rain
  REAL, allocatable :: prcrate1(:,:,:)      ! 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 :: radfrc1(:,:,:)       ! Radiation forcing (K/s)
  REAL, allocatable :: radsw1 (:,:)         ! Solar radiation reaching the surface
  REAL, allocatable :: rnflx1 (:,:)         ! Net radiation flux absorbed by surface

  REAL, allocatable :: usflx1 (:,:)         ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, allocatable :: vsflx1 (:,:)         ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, allocatable :: ptsflx1(:,:)         ! Surface heat flux (K*kg/(m*s**2))
  REAL, allocatable :: qvsflx1(:,:)         ! Surface moisture flux (kg/(m**2*s))

  REAL, allocatable :: ubar1 (:,:,:)        ! Base state u-velocity (m/s).
  REAL, allocatable :: vbar1 (:,:,:)        ! Base state v-velocity (m/s).
  REAL, allocatable :: ptbar1(:,:,:)        ! Base state potential temperature (K)
  REAL, allocatable :: pbar1 (:,:,:)        ! Base state pressure (Pascal).
  REAL, allocatable :: rhobar1(:,:,:)       ! Base state air density (kg/m**3)
  REAL, allocatable :: qvbar1(:,:,:)        ! Base state water vapor specific humidity
                                            ! (kg/kg).

  REAL, allocatable :: x1    (:)            ! The x-coord. of the physical and
                                            ! computational grid. Defined at u-point.
  REAL, allocatable :: y1    (:)            ! The y-coord. of the physical and
                                            ! computational grid. Defined at v-point.
  REAL, allocatable :: z1    (:)            ! The z-coord. of the computational grid.
                                            ! Defined at w-point on the staggered grid.
  !wdt update
  REAL, allocatable :: x1_out(:)
  REAL, allocatable :: y1_out(:)
  REAL, allocatable :: zp1   (:,:,:)        ! The physical height coordinate defined at
                                            ! w-point on the staggered grid.
  REAL, allocatable :: hterain1(:,:)        ! Terrain height (m)
!wdt update
  REAL, allocatable :: htrn1orig(:,:)        ! Terrain height (m)
  REAL, allocatable :: j11   (:,:,:)        ! Coordinate transformation Jacobian -d(zp)/dx.
  REAL, allocatable :: j21   (:,:,:)        ! Coordinate transformation Jacobian -d(zp)/dy.
  REAL, allocatable :: j31   (:,:,:)        ! Coordinate transformation Jacobian  d(zp)/dz.

  REAL, allocatable :: tem11 (:,:,:)        ! Work array
  REAL, allocatable :: tem21 (:,:,:)        ! Work array

  REAL, allocatable :: wgtsx(:,:)  ! Weight for interpolation in x-dir for scalar points
  REAL, allocatable :: wgtsy(:,:)  ! Weight for interpolation in y-dir for scalar points
  REAL, allocatable :: wgtux(:,:)  ! Weight for interpolation in x-dir for u points
  REAL, allocatable :: wgtvy(:,:)  ! Weight for interpolation in y-dir for v points

  REAL, allocatable :: wgtz(:,:,:,:)  ! Weight for interpolation in z-dir for v points

  INTEGER, allocatable :: isx(:),jsy(:),iux(:),jvy(:),kz(:,:,:)

  REAL, allocatable :: zp1d1 (:)       ! Temporary array
  REAL, allocatable :: dzp1d1(:)       ! Temporary array
!
!-----------------------------------------------------------------------
!
!  More arrays.
!
!-----------------------------------------------------------------------
!
  REAL, allocatable :: xs (:), ys (:) ! x,y coord for scalar points
  REAL, allocatable :: xs1(:), ys1(:) ! x,y coord for scalar points

  REAL, allocatable :: xs_2d(:,:), ys_2d(:,:) ! used by subroutine extmnsnd

  REAL, allocatable :: temx1yz  (:,:,:)
  REAL, allocatable :: temx1y1z (:,:,:)
  REAL, allocatable :: temx1y1zb(:,:,:)

  REAL, allocatable :: temz1d1(:),temz1d2(:),temz1d3(:),temz1d4(:),     &
        temz1d5(:),temz1d6(:),temz1d7(:)

  REAL, allocatable :: zsnd(:),ptsnd(:),qvsnd(:),psnd(:),               &
        ubsnd(:),vbsnd(:)

  REAL, allocatable :: ptpsfc(:,:),ppsfc (:,:),qvsfc(:,:),qcsfc(:,:)
  REAL, allocatable :: qrsfc (:,:),qisfc (:,:),qssfc(:,:),qhsfc(:,:)
  REAL, allocatable :: tkesfc(:,:),kmhsfc(:,:),kmvsfc(:,:)
  REAL, allocatable :: ptbsfc(:,:),pbsfc(:,:),qvbsfc(:,:)
  REAL, allocatable :: htrnx1y1(:,:),radsfc(:,:)
  INTEGER, allocatable :: ktrnx1y1(:,:)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
  INTEGER, allocatable :: ix(:,:),jy(:,:)
  REAL, allocatable :: xw(:,:),yw(:,:)

!
!-----------------------------------------------------------------------
!  nx, ny, nz: Dimensions of input grid.
!  nx1, ny1, nz1: Dimensions of output grid.
!-----------------------------------------------------------------------
!
  INTEGER :: nx       ! Number of input grid points in the x-direction
  INTEGER :: ny       ! Number of input grid points in the y-direction
  INTEGER :: nz       ! Number of input grid points in the z-direction

  INTEGER :: nx1      ! Number of output grid points in the x-direction
  INTEGER :: ny1      ! Number of output grid points in the y-direction
  INTEGER :: nz1      ! Number of output grid points in the z-direction

  INTEGER :: nxyz,nxy,nxyz1,nxy1
!
!-----------------------------------------------------------------------
!
!  Soil types.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nstypin           ! Number of soil types in input data
  INTEGER :: nstyp1            ! Number of soil types for ouput data
!
!-----------------------------------------------------------------------
!
!  Misc
!
!-----------------------------------------------------------------------
!
  INTEGER :: lvlprof           ! Number of levels in 1-d average sounding
!
!-----------------------------------------------------------------------
!
  INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
  PARAMETER (nhisfile_max=200)
  CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
  INTEGER :: ireturn

  NAMELIST /output_dims/ nx1,ny1,nz1

  INTEGER :: max_plevels
  PARAMETER (max_plevels=1000)
  REAL :: plevels(max_plevels)
                       ! Array contain the values of pressure levels to
                       ! which fields will be interpolated.
!
!-----------------------------------------------------------------------
!
!  Parameters for the output grid.
!
!  Note:
!
!  Given nx1, ny1 and nz1, the physical domain size of the refined
!  grid will be xl1=(nx1-3)*dx1 by yl1=(ny1-3)*dy1 by zh1=(nz1-3)*dz1.
!  Dx1, dy1 and dz1 are the grid intervals of the refined grid.
!
!-----------------------------------------------------------------------
!
  REAL :: xorig1, yorig1, zorig1
  REAL :: xctr1 , yctr1
  REAL :: ctrlat1, ctrlon1
  REAL :: origlat, origlon

  REAL :: dx1,dy1     ! Grid intervals of the refined grid.

  INTEGER :: strhopt1 ! Vertical grid stretching option.
                      ! = 0, no stretching in vertical.
                      ! >= 1, with stretching in vertical.
  REAL :: dz1         ! Average grid spacing in vertical direction in
                      ! transformed computational space (m).
  REAL :: dzmin1      ! Minimun grid spacing in vertical direction in
                      ! physcal space (m).

  REAL :: zrefsfc1    ! The reference height of the surface
                      ! (ground level) (m)

  REAL :: dlayer11    ! The depth of the lower layer with uniform
                      ! (dz=dzmin) vertical spacing (m)

  REAL :: dlayer21    ! The depth of the mid layer with stetched
                      ! vertical spacing (m)

  REAL :: strhtune1   ! Tuning parameter for stretching option 2
                      ! A Value between 0.2 and 5.0 is appropriate.
                      ! A larger value gives a more linear stretching.

  REAL :: zflat1      ! The height at which the grid levels
                      ! becomes flat in the terrain-following
                      ! coordinate transformation (m).
!
!-----------------------------------------------------------------------
!
!  Terrain option parameters for the new grid:
!
!-----------------------------------------------------------------------
!

  INTEGER :: ternopt1        ! Model terrain option.
                             ! = 0, no terrain, the ground is flat;
                             ! = 1, bell-shaped mountain;
                             ! = 2, terrain data read from terrain data
                             !      base (not implemented yet)
!wdt-williams 2002-01-17 GMB ntmerge
  INTEGER :: ntmerge         ! Number of zones over which to merge original 
                             ! and new terrain at the bondaries 
                             ! (used when ternopt1=4).
  INTEGER :: ternfmt1        ! Terrain data file format.
  INTEGER :: mntopt1         ! Option for choosing idealized mountain
                             ! type.
  REAL :: hmount1         ! The mountain height (m)
  REAL :: mntwidx1        ! The half-width of bell-shaped mountain
                          ! in x-dir.
  REAL :: mntwidy1        ! The half-width of bell-shaped mountain
                          ! in y-dir.
  REAL :: mntctrx1        ! The x coordinate of the bell-shaped
                          ! mountain center.
  REAL :: mntctry1        ! The y coordinate of the bell-shaped
                          ! mountain center.
  CHARACTER (LEN=132) :: terndta1  ! Name of the terrain data file

  REAL :: zflat11,za,zb
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: nunit
  INTEGER :: i, j, k
  REAL :: amin, amax
  REAL :: tmp1,tmp2
  REAL :: radnd,pi2,z0
  INTEGER :: k0,kk
  CHARACTER (LEN=132) :: basdmpfn
  INTEGER :: lbasdmpf
  CHARACTER (LEN=132) :: ternfn,sfcoutfl,soiloutfl,temchar
  INTEGER :: lternfn,lfn
  INTEGER :: iss

!wdt-williams 2002-01-17 GMB ntmerge
  DOUBLE PRECISION :: ntmergeinv, mfac
  INTEGER :: idist
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'indtflg.inc'
  INCLUDE 'alloc.inc'
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
  INCLUDE 'exbc.inc'

  INTEGER :: houtfmt, nchin, nchout
  CHARACTER (LEN=132) :: filename

  INTEGER :: grdbas

  REAL :: time
  !wdt update
  INTEGER :: vroutcnt
  DATA    vroutcnt /0/

  INTEGER :: nfile, length, lenstr
  CHARACTER (LEN=132) :: timsnd
  CHARACTER (LEN=80) :: new_runname
  INTEGER :: tmstrln
  REAL :: xeps, yeps
  REAL :: ctrx,ctry,swx,swy,alatpro(2),sclf,dxscl,dyscl
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
  INTEGER :: bglopt,ntagopt
  REAL :: misvalue,aghght

  REAL :: snddelz,delz,dtdz,tbark1,tbark,ttotk1,ttotk,lnpbar,lnptot

  INTEGER :: i1mn,i1mx,j1mn,j1mx

  INTEGER :: npoint_below_ground ! flag indicating if any point in the
                                 ! output grid terrain is below input grid terrain
  INTEGER :: redo_base_state
  REAL :: zpmin,zpmax,zpmin1,pbartop, tem,tema,a,b,c,d

!
!-----------------------------------------------------------------------
!
!  namelist Declarations:
!
!-----------------------------------------------------------------------
!
  NAMELIST /jobname/ runname

  INTEGER :: z_or_p   ! Control parameter denoting if the output grid is
                      ! in height or pressure coordinates.
  INTEGER :: xy_or_ll ! Control parameter denoting if (xctr1,yctr1) or
                      ! (ctrlat1,ctrlon1) are to be used to specify the new
                      ! grid center.
  INTEGER :: intrphopt  ! Option for horizontal interpolation
                        ! = 1 linear, =2 quadratic
  INTEGER :: intrpvopt  ! Option for vertical interpolation
                        ! = 1 linear, =2 quadratic

  INTEGER :: snap_to_grid, same_res
  INTEGER :: istatus, ib,jb

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
  INTEGER, PARAMETER :: MAX_GRD=200
  INTEGER :: noutgrds
  REAL :: xctr_grd(MAX_GRD), yctr_grd(MAX_GRD)
  REAL :: clat_grd(MAX_GRD), clon_grd(MAX_GRD)
  CHARACTER (LEN=80 ) :: name_grd(MAX_GRD)
  INTEGER :: ng

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
  REAL ::  ctrlat_sv, ctrlon_sv, latitud_sv, dx_sv, dy_sv, dz_sv
  REAL :: dzmin_sv, zrefsfc_sv, dlayer1_sv, dlayer2_sv, zflat_sv, strhtune_sv
  INTEGER :: strhopt_sv, nstyp_sv

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
  NAMELIST /newgrid/ z_or_p,xy_or_ll,strhopt1,xctr1,yctr1,              &
          ctrlat1,ctrlon1,snap_to_grid,same_res,dx1,dy1,dz1,dzmin1,     &
          zrefsfc1,dlayer11,dlayer21,strhtune1,zflat1,plevels,nstyp1,   &
          noutgrds,xctr_grd,yctr_grd,clat_grd,clon_grd,name_grd

!wdt-williams 2002-01-17 GMB ntmerge
  NAMELIST /newterrain/ ternopt1,mntopt1,hmount1,mntwidx1,mntwidy1,     &
            mntctrx1,mntctry1,terndta1,ternfmt1,ntmerge

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
  NAMELIST /bgloption/ bglopt,misvalue,intrphopt,intrpvopt,ntagopt,aghght

  NAMELIST /sfc_data/ sfcdat,styp,vtyp,lai0,roufns0,veg0,               &
            sfcdtfl,sfcfmt

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
  NAMELIST /output/ dirname,exbcdmp,hdmpfmt,grbpkbit,hdfcompr,          &
            grdout,basout,varout,mstout,rainout,prcout,iceout,          &
            tkeout, trbout,sfcout,landout,radout,flxout,                &
            qcexout,qrexout,qiexout,qsexout,qhexout,                    &
            totout,filcmprs,readyfl,sfcdmp,soildmp,terndmp,ngbrz,zbrdmp
!wdt update
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  REAL :: cpu0, cpu_time, f_cputime
  REAL :: xctr1_old,yctr1_old,xsw,ysw,xsw1,ysw1,xsw1_old,ysw1_old
  REAL :: ctrlat1_old, ctrlon1_old

!wdt update
  REAL :: rhmax = 1.0
  REAL :: qvsat

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
  REAL :: hght0, dhght, fac

!wdt update
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

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

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  cpu0 = f_cputime()

  WRITE(6,'(/9(/2x,a)/)')                                               &
     '###############################################################', &
     '###############################################################', &
     '###                                                         ###', &
     '###                Welcome to ARPSINTRP                     ###', &
     '###      This program converts the history dump data        ###', &
     '###      sets generated by ARPS, between various formats.   ###', &
     '###                                                         ###', &
     '###############################################################', &
     '###############################################################'

!
!-----------------------------------------------------------------------
!
!  Read in the input parameters.
!
!-----------------------------------------------------------------------
!
  CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)

  lengbf = len_trim(grdbasfn)

  CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf),                    &
       nx,ny,nz,nstypin, ireturn)
  nstypin = max(1,nstypin)

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

  WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz

  READ(5,output_dims, END=105)

  WRITE(6,'(/a,a/)')                                                    &
      'NAMELIST block intrp_dims sucessfully read.'

  WRITE(6,'(3(a,i5))') 'nx1=',nx1,', ny1=',ny1,', nz1=',nz1

  GO TO 10

  105   WRITE(6,'(/a,a/)')                                              &
          'Error reading NAMELIST block intrp_dims. ',                  &
          'Program ARPSINTRP stopped.'
  STOP 1

  10    CONTINUE

  WRITE (6,'(/a,g13.3/)') "Current memory allocation (in words):",      &
               current_memory_use
!

  IF (max_plevels < nz1) THEN
    WRITE (6,*) "ARPSINTRP: ERROR, nz1 < max_plevels.  Stopping."
    WRITE (6,*) "nz1 =",nz1,"  max_plevels = ",max_plevels
    STOP 1
  END IF

  mgrid = 1
  nestgrd = 0
!
!-----------------------------------------------------------------------
!
!  Set the default parameters
!
!-----------------------------------------------------------------------
!
  runname   = 'intrp'

  z_or_p    = 1
  xy_or_ll  = 1
  xctr1     = 512000.0
  yctr1     = 512000.0
  ctrlat1   = 34.0
  ctrlon1   = -98.0
  dx1       =  4000.0
  dy1       =  4000.0
  dz1       =  350.0
  strhopt1  = 2
  dzmin1    = 20.0
  zrefsfc1  = 0.0
  dlayer11  = 0.0
  dlayer21  = 1.0E5
  strhtune1 = 1.0
  zflat1    = 1.0E5
  nstyp1    = -1

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
  noutgrds = 0
  xctr_grd = 0
  yctr_grd = 0
  clat_grd = ctrlat1
  clon_grd = ctrlon1
  name_grd = 'NULL'

  ternopt1  = 3
!wdt-williams 2002-01-17 GMB ntmerge
  ntmerge = 1
  mntopt1   = 1
  hmount1   =     0.000
  mntwidx1  = 10000.000
  mntwidy1  = 10000.000
  mntctrx1  = 10000.000
  mntctry1  = 10000.000
  terndta1  = 'arpstern.dat'
  ternfmt1  = 0

  bglopt    = 1
  misvalue  = -9999.0
  intrphopt = 1
  intrpvopt = 1
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
  ntagopt = 0
  aghght = 500.0

  sfcdat  = 1
  styp    = 3
  vtyp    = 10
  lai0    = 0.31
  roufns0 = 0.1
  veg0    = 0.0
  sfcdtfl = 'arpssfc.data'
  sfcfmt  = 0

  dirname   = './'
  exbcdmp   = 1
  hdmpfmt   = 1
  grbpkbit  = 16
  hdfcompr  = 0
  filcmprs  = 0
  readyfl   = 1
  basout    = 0
  grdout    = 0
  varout    = 1
  mstout    = 1
  iceout    = 1
  tkeout    = 1
  trbout    = 0
  rainout   = 0
  sfcout    = 0
  snowout   = 0
  landout   = 0
  qcexout   = 0
  qrexout   = 0
  qiexout   = 0
  qsexout   = 0
  qhexout   = 0

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 
  ngbrz = 5
  zbrdmp = 10000.0

  sfcdmp    = 1
  soildmp   = 1
  terndmp   = 1
  same_res  = 0
  snap_to_grid = 0

  READ (5,jobname,END=100)
  WRITE(6,'(/a/)') 'Sucessfully read namelist block JOBNAME.'

  WRITE(6,'(/2x,a,a)') 'The name of this run is: ', runname
  new_runname = runname
!
!-----------------------------------------------------------------------
!
!  Set the output grid and the variable control parameters
!
!-----------------------------------------------------------------------
!
  READ (5,newgrid)
  WRITE(6,'(/a/)') 'Sucessfully read namelist block NEWGRID.'

  PRINT*
  PRINT*,' Input parameters for the new refined grid:'
  PRINT*

  PRINT*,' z_or_p= ', z_or_p
  PRINT*,' xy_or_ll= ', xy_or_ll
  PRINT*,' Input dx1:'
  PRINT*,' dx1    = ',dx1
  PRINT*,' dy1    =  ',dy1
  PRINT*,' strhopt1= ',strhopt1
  PRINT*,' dz1    =  ',dz1
  PRINT*,' dzmin1 =  ',dzmin1
  PRINT*,' xctr1  =  ',xctr1
  PRINT*,' yctr1  =  ',yctr1
  PRINT*,' ctrlat1  =  ',ctrlat1
  PRINT*,' ctrlon1  =  ',ctrlon1
  PRINT*,'plevels =', (plevels(k),k=1,nz1)
  PRINT*,' nstyp1   =  ',nstyp1

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
  PRINT*,'noutgrds = ',noutgrds
  IF (noutgrds > MAX_GRD) THEN
    WRITE (*,*) "ERROR, noutgrds greater than MAX_GRD. ",  &
       "Set noutgrds to less than or equal to ",MAX_GRD," ABORTING!"
    STOP
  ENDIF
  IF (noutgrds > 0) THEN
    DO ng = 1,noutgrds
      PRINT*,'xctr_grd(',ng,') = ',xctr_grd(ng)
      PRINT*,'yctr_grd(',ng,') = ',yctr_grd(ng)
      PRINT*,'clat_grd(',ng,') = ',clat_grd(ng)
      PRINT*,'clon_grd(',ng,') = ',clon_grd(ng)
      PRINT*,'name_grd(',ng,') = ',name_grd(ng)
    END DO
  ENDIF

!
!-----------------------------------------------------------------------
!
!  Specify the terrain initialization option parameters:
!
!-----------------------------------------------------------------------
!

  READ (5,newterrain,END=100)
  WRITE(6,'(/a)') 'Sucessfully read namelist block NEWTERRAIN.'

  WRITE(6,'(1x,a,i4)') 'ternopt1 =', ternopt1
!wdt-williams 2002-01-17 GMB ntmerge
  WRITE(6,'(1x,a,i4)') 'ntmerge  =', ntmerge
  WRITE(6,'(1x,a,i4)') 'mntopt1  =', mntopt1
  WRITE(6,'(1x,a,f10.3)') 'hmount1 =',hmount1
  WRITE(6,'(1x,a,f10.3)') 'mntwidx1=',mntwidx
  WRITE(6,'(1x,a,f10.3)') 'mntwidy1=',mntwidy
  WRITE(6,'(1x,a,f10.3)') 'mntctrx1=',mntctrx
  WRITE(6,'(1x,a,f10.3)') 'mntctry1=',mntctry

  lenstr = 132
  CALL strlnth( terndta1, lenstr)
  WRITE(6,'(2x,a,a/)') 'terndta1 = ',terndta1(1:lenstr)
  WRITE(6,'(1x,a,i4)') 'ternfmt1 =', ternfmt1
  ternfmt = ternfmt1


  READ(5,bgloption,END=100)
  WRITE(6,'(/a)') 'Sucessfully read namelist block BLGOPTION.'

  WRITE(6,'(1x,a,i4)')    'bglopt  =', bglopt
  WRITE(6,'(1x,a,f10.3)') 'misvalue=', misvalue
  WRITE(6,'(1x,a,i4)')    'intrphopt  =', intrphopt
  WRITE(6,'(1x,a,i4)')    'intrpvopt  =', intrpvopt
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
  WRITE(6,'(1x,a,i4)')    'ntagopt  =', ntagopt
  WRITE(6,'(1x,a,f10.3)') 'aghght  =', aghght

!
!-----------------------------------------------------------------------
!
!  Set the control parameters for soil data:
!
!-----------------------------------------------------------------------
!
  READ (5,sfc_data,END=100)
  WRITE(6,'(/a/)') 'Sucessfully read namelist block SOIL_DATA.'

  IF (sfcdat .eq. 1) THEN
    nstyp = 1
    nstypin = nstyp
  ELSE
    nstyp = nstypin
  ENDIF
  IF (nstyp1 .eq. -1) nstyp1 = nstyp

  WRITE (6,'(a,i3)') "Number of soil types to be used for input data:",nstypin
!
!-----------------------------------------------------------------------
!
!  Set the control parameters for output:
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)')                                                        &
      ' Reading in control parameters for the output data files..'

  READ (5,output,END=100)
  WRITE(6,'(a)') 'Sucessfully read namelist block OUTPUT.'

  houtfmt = hdmpfmt

!  IF( houtfmt.eq.0 ) THEN
!    write(6,'(/1x,a/)') 'Output format is 0, no data is dumped.'
!    STOP 9103
!  ENDIF

  IF ( houtfmt == 10 .AND. grbpkbit <= 0 ) THEN
    WRITE(6,'(a,a,i2/a)')                                               &
        'The bit width for packing GRIB data was invalid, ',            &
        'The old value was ', grbpkbit, 'Reset it to 16 bits'
    grbpkbit = 16
  END IF

  totout = 1

  IF (sfcout == 1) THEN
    snowout = 1
  END IF

  GO TO 15

  100   WRITE(6,'(a)')                                                  &
          'Error reading NAMELIST file. Program ARPSINTRP stopped.'
  STOP 9104

  15    CONTINUE

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

  lengbf=len_trim(grdbasfn)
  WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)
!
!-----------------------------------------------------------------------
!
!  Allocate arrays.
!
!-----------------------------------------------------------------------
!
  WRITE (6,'(/a,g13.3/)') "Current memory allocation (in words):",      &
               current_memory_use

!
  lvlprof=MAX(601,nz)
  nxyz = nx*ny*nz
  nxy = nx*ny

  allocate(  CALL alloc_status_accounting(istatus,nxyz,'u')
  allocate(  CALL alloc_status_accounting(istatus,nxyz,'v')
  allocate(  CALL alloc_status_accounting(istatus,nxyz,'w')
  allocate(ptprt(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'ptprt')
  allocate(pprt(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'pprt')

  allocate(qv(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qv')
  allocate(qc(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qc')
  allocate(qr(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qr')
  allocate(qi(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qi')
  allocate(qs(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qs')
  allocate(qh(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qh')
  allocate(tke(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'tke')
  allocate(kmh(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'kmh')
  allocate(kmv(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'kmh')

  allocate(soiltyp(nx,ny,nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*nstypin,'soiltyp')
  allocate(stypfrct(nx,ny,nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*nstypin,'stypfrct')
  allocate(vegtyp(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'vegtyp ')
  allocate(lai(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'lai    ')
  allocate(roufns(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'roufns ')
  allocate(veg(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'veg    ')
  allocate(ndvi(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'ndvi   ')

  allocate(tsfc(nx,ny,0:nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'tsfc   ')
  allocate(tsoil(nx,ny,0:nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'tsoil  ')
  allocate(wetsfc(nx,ny,0:nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetsfc ')
  allocate(wetdp(nx,ny,0:nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetdp  ')
  allocate(wetcanp(nx,ny,0:nstypin),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*(nstypin+1),'wetcanp')
  allocate(snowdpth(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'snowdpth')

  allocate(raing(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'raing')
  allocate(rainc(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'rainc')
  allocate(prcrate(nx,ny,4),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy*4,'prcrate')

  allocate(radfrc(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'radfrc')
  allocate(radsw(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'radsw ')
  allocate(rnflx(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'rnflx ')
  allocate(usflx(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'usflx ')
  allocate(vsflx(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'vsflx ')
  allocate(ptsflx(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'ptsflx')
  allocate(qvsflx(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'qvsflx')

  allocate(ubar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'ubar')
  allocate(vbar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'vbar')
  allocate(wbar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'wbar')
  allocate(ptbar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'ptbar')
  allocate(pbar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'pbar')
  allocate(rhobar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'rhobar')
  allocate(qvbar(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qvbar')

  allocate(  CALL alloc_status_accounting(istatus,nx,'x')
  allocate(  CALL alloc_status_accounting(istatus,ny,'y')
  allocate(  CALL alloc_status_accounting(istatus,nz,'z')
  allocate(zp(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'zp')
  allocate(hterain(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'hterain')

  allocate(uprt(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'uprt')
  allocate(vprt(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'vprt')
  allocate(qvprt(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'qvprt')

  allocate(tem1(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'tem1')
  allocate(tem2(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'tem2')
  allocate(tem3(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'tem3')
  !wdt update
  allocate(tem4(nx,ny,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'tem4')

  allocate(u1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'u1')
  allocate(v1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'v1')
  allocate(w1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'w1')
  allocate(ptprt1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'ptprt1')
  allocate(pprt1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz,'pprt1')

  nxyz1 = nx1*ny1*nz1
  nxy1  = nx1*ny1

  print*,'Start allocating arrays for output grid.'

  allocate(qv1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qv1')
  qv1 = 0.0
  allocate(qc1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qc1')
  qc1 = 0.0
  allocate(qr1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qr1')
  qr1 = 0.0
  allocate(qi1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qi1')
  qi1 = 0.0
  allocate(qs1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qs1')
  qs1 = 0.0
  allocate(qh1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qh1')
  qh1 = 0.0
  allocate(tke1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'tke1')
  tke1= 0.0
  allocate(kmh1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'kmh1')
  kmh1= 0.0
  allocate(kmv1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'kmv1')
  kmv1= 0.0

  allocate(soiltyp1(nx1,ny1,nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*nstyp1,'soiltyp1')
  allocate(stypfrct1(nx1,ny1,nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*nstyp1,'stypfrct1')
  soiltyp1 = 0
  stypfrct1 = 0.0
  allocate(vegtyp1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'vegtyp1')
  vegtyp1 = 0
  allocate(lai1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'lai1')
  lai1 = 0.0
  allocate(roufns1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'roufns1')
  roufns1 = 0.0
  allocate(veg1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'veg1')
  veg1 = 0.0

  allocate(tsfc1(nx1,ny1,0:nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'tsfc1')
  allocate(tsoil1(nx1,ny1,0:nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'tsoil1')
  allocate(wetsfc1(nx1,ny1,0:nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetsfc1')
  allocate(wetdp1(nx1,ny1,0:nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetdp1')
  allocate(wetcanp1(nx1,ny1,0:nstyp1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'wetcanp1')
  tsfc1 = 0.0
  tsoil1 = 0.0
  wetsfc1 = 0.0
  wetdp1 = 0.0
  wetcanp1 = 0.0
  allocate(snowdpth1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*(nstyp1+1),'snowdpth1')
  snowdpth1 = 0.0

  allocate(raing1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'raing1')
  raing1 = 0.0
  allocate(rainc1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'rainc1')
  rainc1 = 0.0
  allocate(prcrate1(nx1,ny1,4),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1*4,'prcrate1')
  prcrate1 = 0.0

  allocate(radfrc1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'radfrc1')
  radfrc1 = 0.0
  allocate(radsw1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'radsw1')
  radsw1 = 0.0
  allocate(rnflx1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'rnflx1')
  rnflx1 = 0.0

  allocate(usflx1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'usflx1')
  usflx1 = 0.0
  allocate(vsflx1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'vsflx1')
  vsflx1 = 0.0
  allocate(ptsflx1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'ptsflx1')
  ptsflx1 = 0.0
  allocate(qvsflx1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qvsflx1')
  qvsflx1 = 0.0

  allocate(ubar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'ubar1')
  ubar1 = 0.0
  allocate(vbar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'vbar1')
  vbar1 = 0.0
  allocate(ptbar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'pbar1')
  ptbar1 = 0.0
  allocate(pbar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'ptbar1')
  pbar1 = 0.0
  allocate(rhobar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'rhobar1')
  rhobar1 = 0.0
  allocate(qvbar1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'qvbar1')
  qvbar1 = 0.0

  allocate(x1(nx1),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1,'x1')
  x1 = 0.0
  allocate(y1(ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1,'y1')
  y1 = 0.0
  !wdt update
  allocate(x1_out(nx1),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1,'x1_out')
  x1_out = 0.0
  allocate(y1_out(ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1,'y1_out')
  y1_out = 0.0
  allocate(z1(nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nz1,'z1')
  z1 = 0.0
  allocate(zp1(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'zp1')
  zp1 = 0.0
  allocate(hterain1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'hterain1')
  hterain1 = 0.0
!wdt update
  allocate(htrn1orig(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'htrn1orig')
  htrn1orig = 0.0
  allocate(j11(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'j11')
  j11 = 0.0
  allocate(j21(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'j21')
  j21 = 0.0
  allocate(j31(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'j31')
  j31 = 1.0

  allocate(tem11(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'tem11')
  tem11 = 0.0
  allocate(tem21(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'tem21')
  tem21 = 0.0

  allocate(wgtsx(nx1,3),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1*3,'wgtsx')
  wgtsx = 0.0
  allocate(wgtsy(ny1,3),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1*3,'wgtsy')
  wgtsy = 0.0
  allocate(wgtux(nx1,3),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1*3,'wgtux')
  wgtux = 0.0
  allocate(wgtvy(ny1,3),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1*3,'wgtvy')
  wgtvy = 0.0
  allocate(wgtz(nx1,ny1,nz1,3),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1*3,'wgtz')
  wgtz  = 0.0

  allocate(isx(nx1),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1,'isx')
  isx = 0.0
  allocate(jsy(ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1,'jsy')
  jsy = 0.0
  allocate(iux(nx1),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1,'iux')
  iux = 0.0
  allocate(jvy(ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1,'jvy')
  jvy = 0.0
  allocate(kz(nx1,ny1,nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxyz1,'kz')
  kz = 0.0
  allocate(zp1d1(nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nz1,'zp1d1')
  zp1d1 = 0.0
  allocate(dzp1d1(nz1),stat=istatus)
  CALL alloc_status_accounting(istatus,nz1,'dzp1d1')
  dzp1d1 = 0.0

  allocate(xs(nx ),stat=istatus)
  CALL alloc_status_accounting(istatus,nx,'xs')
  xs = 0.0
  allocate(ys(ny ),stat=istatus)
  CALL alloc_status_accounting(istatus,ny,'ys')
  ys = 0.0
  allocate(xs1(nx1),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1,'xs1')
  xs1 = 0.0
  allocate(ys1(ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,ny1,'ys1')
  ys1 = 0.0

  allocate(xs_2d(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'xs_2d')
  xs_2d = 0.0
  allocate(ys_2d(nx,ny),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy,'ys_2d')
  ys_2d = 0.0

  allocate(temx1yz(nx1,ny ,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1*ny*nz,'temx1yz')
  temx1yz = 0.0
  allocate(temx1y1z(nx1,ny1,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1*ny1*nz,'temx1y1z')
  temx1y1z = 0.0
  allocate(temx1y1zb(nx1,ny1,nz),stat=istatus)
  CALL alloc_status_accounting(istatus,nx1*ny1*nz,'temx1y1zb')
  temx1y1zb = 0.0

  allocate(temz1d1(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d1')
  temz1d1 = 0.0
  allocate(temz1d2(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d2')
  temz1d2 = 0.0
  allocate(temz1d3(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d3')
  temz1d3 = 0.0
  allocate(temz1d4(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d4')
  temz1d4 = 0.0
  allocate(temz1d5(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d5')
  temz1d5 = 0.0
  allocate(temz1d6(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d6')
  temz1d6 = 0.0
  allocate(temz1d7(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'temz1d7')
  temz1d7 = 0.0

  allocate(zsnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'zsnd')
  zsnd = 0.0
  allocate(ptsnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'ptsnd')
  ptsnd = 0.0
  allocate(qvsnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'qvsnd')
  qvsnd = 0.0
  allocate(psnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'psnd')
  psnd = 0.0
  allocate(ubsnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'ubsnd')
  ubsnd = 0.0
  allocate(vbsnd(lvlprof),stat=istatus)
  CALL alloc_status_accounting(istatus,lvlprof,'vbsnd')
  vbsnd = 0.0

  allocate(ptpsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'ptpsfc')
  ptpsfc = 0.0
  allocate(ppsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'ppsfc')
  ppsfc = 0.0
  allocate(qvsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qvsfc')
  qvsfc = 0.0
  allocate(qcsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qcsfc')
  qcsfc = 0.0
  allocate(qrsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qrsfc')
  qrsfc = 0.0
  allocate(qisfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qisfc')
  qisfc = 0.0
  allocate(qssfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qssfc')
  qssfc = 0.0
  allocate(qhsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qhsfc')
  qhsfc = 0.0
  allocate(tkesfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'tkesfc')
  tkesfc = 0.0
  allocate(kmhsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'kmhsfc')
  kmhsfc = 0.0
  allocate(kmvsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'kmvsfc')
  kmvsfc = 0.0
  allocate(ptbsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'ptbsfc')
  ptbsfc = 0.0
  allocate(pbsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'pbsfc')
  pbsfc = 0.0
  allocate(qvbsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'qvbsfc')
  qvbsfc = 0.0
  allocate(htrnx1y1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'htrnx1y1')
  htrnx1y1 = 0.0
  allocate(radsfc(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'radsfc')
  radsfc = 0.0
  allocate(ktrnx1y1(nx1,ny1),stat=istatus)
  CALL alloc_status_accounting(istatus,nxy1,'ktrnx1y1')
  ktrnx1y1 = 0.0

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
  allocate (ix(nx1,ny1),stat=istatus)
  allocate (jy(nx1,ny1),stat=istatus)
  allocate (xw(nx1,ny1),stat=istatus)
  allocate (yw(nx1,ny1),stat=istatus)

  nxyz = nx*ny*nz

  CALL flzero(u     , nxyz)
  CALL flzero(v     , nxyz)
  CALL flzero(uprt  , nxyz)
  CALL flzero(vprt  , nxyz)
  CALL flzero(w     , nxyz)
  CALL flzero(ptprt , nxyz)
  CALL flzero(pprt  , nxyz)
  CALL flzero(qv    , nxyz)
  CALL flzero(qvprt , nxyz)
  CALL flzero(qc    , nxyz)
  CALL flzero(qr    , nxyz)
  CALL flzero(qi    , nxyz)
  CALL flzero(qs    , nxyz)
  CALL flzero(qh    , nxyz)
  CALL flzero(tke   , nxyz)
  CALL flzero(kmh   , nxyz)
  CALL flzero(kmv   , nxyz)

  CALL flzero(ubar  , nxyz)
  CALL flzero(vbar  , nxyz)
  CALL flzero(wbar  , nxyz)
  CALL flzero(ptbar , nxyz)
  CALL flzero(pbar  , nxyz)
  CALL flzero(rhobar, nxyz)
  CALL flzero(qvbar , nxyz)

  CALL flzero(x, nx)
  CALL flzero(y, ny)
  CALL flzero(z, nz)
  CALL flzero(zp    , nxyz)

  DO iss=0,nstyp
    DO j=1,ny
      DO i=1,nx
        tsfc   (i,j,iss) = 0.0
        tsoil  (i,j,iss) = 0.0
        wetsfc (i,j,iss) = 0.0
        wetdp  (i,j,iss) = 0.0
        wetcanp(i,j,iss) = 0.0
      END DO
    END DO
  END DO

  DO iss=1,nstyp
    DO j=1,ny
      DO i=1,nx
        soiltyp (i,j,iss) = 0
        stypfrct(i,j,iss) = 0.0
      END DO
    END DO
  END DO

  vegtyp  = 0
  snowdpth= 0.0
  lai     = 0.0
  roufns  = 0.0
  veg     = 0.0
  ndvi    = 0.0
  raing   = 0.0
  rainc   = 0.0
  prcrate = 0.0
  prcrate = 0.0
  prcrate = 0.0
  prcrate = 0.0
  radsw   = 0.0
  rnflx   = 0.0
  usflx   = 0.0
  vsflx   = 0.0
  ptsflx  = 0.0
  qvsflx  = 0.0
 
  radfrc = 0.0

!
!-----------------------------------------------------------------------
!
!  Run interpolation.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Loop over data files
!
!-----------------------------------------------------------------------
!
  ireturn = 0

  DO nfile = 1,nhisfile

    filename = hisfile(nfile)

    lenfil=len_trim(filename)
    WRITE(6,'(/a,a,a)')                                                 &
        ' Data set ', filename(1:lenfil) ,' to be processed.'
!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
    102   CONTINUE

    CALL dtaread(nx,ny,nz,nstyp,                                       &
                 hinfmt, nchin,grdbasfn(1:lengbf),lengbf,               &
                 filename(1:lenfil),lenfil,time,                        &
                 x,y,z,zp, uprt ,vprt ,w,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)

    !wdt update moved code block to here
    !wdt begin block
    curtim = time

    IF( hinfmt == 9 .AND. ireturn == 2 ) THEN
      WRITE(6,'(/1x,a/)') 'The end of GrADS file was reached.'
      CLOSE ( nchin )
      CALL retunit( nchin )
      GO TO 9001
    END IF

    IF( ireturn /= 0 ) GO TO 9002            ! Read was unsuccessful
!
!-----------------------------------------------------------------------
!
!  Grid and base related calculations that need to be done only once.
!
!-----------------------------------------------------------------------
!
    IF( nfile == 1) THEN ! Need to do this only once.
!
!-----------------------------------------------------------------------
!
!  Set-up land surface property data, if requested.
!
!-----------------------------------------------------------------------
!
      IF ( sfcdat == 0 ) GO TO 110

      IF ( sfcdat == 1 ) THEN

        nstyp = 1

        DO j=1, ny-1
          DO i=1, nx-1
            soiltyp(i,j,1) = styp
            vegtyp (i,j) = vtyp
            lai    (i,j) = lai0
            roufns (i,j) = roufns0
            veg    (i,j) = veg0
          END DO
        END DO

      ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN
!
!-----------------------------------------------------------------------
!
!  Read the surface property data from file sfcdtfl (passed
!  in globcst.inc).
!
!-----------------------------------------------------------------------
!
        CALL readsfcdt( nx,ny,nstyp,sfcdtfl,dx,dy,                     &
             mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,       &
             soiltyp,stypfrct,vegtyp,lai,roufns,veg,ndvi )

      ELSE IF(sfcdat == 3 .AND. landin == 1) THEN

        WRITE(6,'(1x,a/,a/)')                                           &
            'Surface property data in the history file was used.',      &
            'Data in ',sfcdtfl,' ignored.'

      ELSE

        WRITE(6,'(1x,a,i3,a/)')                                         &
            'Invalid surface data input option. sfcdat =',sfcdat,       &
            '. Program stopped in ARPSINTRP'
        STOP

      END IF  ! sfcdat option

      IF ( nstyp == 1 ) THEN
        DO j=1,ny
          DO i=1,nx
            stypfrct(i,j,1) = 1.0
          END DO
        END DO
      END IF

      110     CONTINUE

      DO j=1,ny
        DO i=1,nx
          hterain(i,j) = zp(i,j,2)
        END DO
      END DO

      dx = x(2) - x(1)
      dy = y(2) - y(1)

      xorig = x(2)
      yorig = y(2)
      zorig = z(2)

      WRITE(6,'(1x,a,3f15.3)') 'xorig, yorig, zorig =',                 &
           xorig, yorig, zorig

      ebc = 3
      wbc = 3
      sbc = 3
      nbc = 3

      dx = x(2)-x(1)
      dy = y(2)-y(1)
      dz = z(2)-z(1)

    END IF
    !wdt end block
!
!-----------------------------------------------------------------------
!
!  Set up the map projection of the input grid.
!
!-----------------------------------------------------------------------
!
    alatpro(1)=trulat1
    alatpro(2)=trulat2

    dx = x(2)-x(1)
    dy = y(2)-y(1)
    IF( sclfct /= 1.0) THEN
      sclf  = 1.0/sclfct
      dxscl = dx*sclf
      dyscl = dy*sclf
    ELSE
      sclf  = 1.0
      dxscl = dx
      dyscl = dy
    END IF

    PRINT*,mapproj,sclf,alatpro,trulon,ctrlat,ctrlon

    !wdt update: moved code block
    IF (same_res == 1) THEN
      IF(dx /= dx1 .OR. dy /= dy1) THEN
        WRITE (*,*) "WARNING: resetting dx1 and/or dy1 to match input grid."
      ENDIF
      dx1 = dx
      dy1 = dy
    ENDIF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
    ng = 0
500 CONTINUE ! Loop here for each output grid (noutgrds)    
    IF (noutgrds > 0) THEN
      ng = ng + 1
      ctrlat1 = clat_grd(ng)
      ctrlon1 = clon_grd(ng)
      xctr1 = xctr_grd(ng)
      yctr1 = yctr_grd(ng)
    ENDIF

    CALL setmapr( mapproj,sclf,alatpro,trulon )
    CALL lltoxy( 1,1, ctrlat,ctrlon, ctrx, ctry )
    swx = ctrx - (FLOAT(nx-3)/2.)*dxscl
    swy = ctry - (FLOAT(ny-3)/2.)*dyscl
    CALL setorig( 1, swx, swy) ! set up the model origin to the coord.

    CALL xytoll(1,1,(FLOAT(nx-3)/2.)*dx,(FLOAT(ny-3)/2.)*dy,            &
         tmp1,tmp2)
    PRINT*,'ctrlat=',tmp1,', ctrlon=',tmp2

    IF (xy_or_ll == 1) THEN
      CALL xytoll(1,1,xctr1,yctr1, ctrlat1,ctrlon1)
    ELSE IF (xy_or_ll == 2) THEN
      CALL lltoxy(1,1, ctrlat1,ctrlon1, xctr1,yctr1)
    ELSE
      WRITE (*,*) "ERROR: xy_or_ll =",xy_or_ll," not a valid value."
      STOP
    END IF


    IF (same_res == 1) THEN
      IF(dx /= dx1 .OR. dy /= dy1) THEN
        WRITE (*,*) "WARNING: resetting dx1 and/or dy1 to match input grid."
      ENDIF
      dx1 = dx
      dy1 = dy
    ENDIF

    ctrlat1_old = ctrlat1
    ctrlon1_old = ctrlon1

    IF( snap_to_grid /= 0 .OR. same_res == 1) THEN

      xsw = x(1)
      ysw = y(1)
      xsw1 = xctr1-(nx1-1)*dx1*0.5
      ysw1 = yctr1-(ny1-1)*dy1*0.5

      xsw1_old = xsw1
      ysw1_old = ysw1

      PRINT*,'mod =', MOD(dx,dx1),MOD(dy,dy1)

      IF(((dx >= dx1 .AND.                                              &
            ABS((MOD(dx,dx1)-dx1)*MOD(dx,dx1)) <= 0.001*dx1).OR.        &
            (dx <= dx1 .AND.                                            &
            ABS((MOD(dx1,dx)-dx )*MOD(dx1,dx)) <= 0.001*dx )).AND.      &
            ((dy >= dy1 .AND.                                           &
            ABS((MOD(dy,dy1)-dy1)*MOD(dy,dy1)) <= 0.001*dy1).OR.        &
            (dy <= dy1 .AND.                                            &
            ABS((MOD(dy1,dy)-dy )*MOD(dy1,dy)) <= 0.001*dy ) )) THEN
        PRINT*,'is multiple'

        IF(snap_to_grid == 1) THEN
          xsw1 = xsw + nint( (xsw1_old-xsw)/dx1 )*dx1
          ysw1 = ysw + nint( (ysw1_old-ysw)/dy1 )*dy1
        ELSE IF(snap_to_grid == 2) THEN
          xsw1 = xsw + nint( (xsw1_old-xsw)/dx )*dx
          ysw1 = ysw + nint( (ysw1_old-ysw)/dy )*dy
        END IF

        PRINT*,'xsw,ysw,xsw1,ysw1=',xsw,ysw,xsw1,ysw1
        PRINT*,'xsw1_old,ysw1_old=',xsw1_old,ysw1_old

        IF((snap_to_grid /= 0) .AND. (                                  &
              (ABS((xsw1_old - xsw1)) > 0.00001*dx1).OR.                &
              (ABS((ysw1_old - ysw1)) > 0.00001*dy1) )) THEN
          xctr1_old = xctr1
          yctr1_old = yctr1

          xctr1 = xsw1+(nx1-1)*dx1*0.5
          yctr1 = ysw1+(ny1-1)*dy1*0.5
!
!    Note there is a possibility that the shifted grid is out of bound
!    The checking is done later.
!

          WRITE(6,'(/1x,a/)')                                           &
              '################################################################'
          WRITE(6,'(/1x,a/,2(/1x,a)/)')                                 &
              '                    ATTENTION:',                         &
              'xctr1 and/or yctr1 reset so that some of the output ',   &
              'grid lines match those of input grid.'

          WRITE(6,'(1x,2(a,f17.5))')'Original xctr1_old=',xctr1_old,    &
                                  ', yctr1_old=',yctr1_old
          WRITE(6,'(1x,2(a,f17.5))')'New      xctr1   = ',xctr1,        &
                                  ', yctr1   = ',yctr1
          WRITE(6,'(/1x,a/)')                                           &
              '################################################################'

          CALL xytoll(1,1,xctr1,yctr1, ctrlat1,ctrlon1)
        END IF
      ELSE
        same_res = 0  ! ???
      END IF

    END IF

    IF( nz1 /= nz .AND. same_res == 1) THEN
      WRITE (*,*) "WARNING: nz1 (",nz1,") does not match nz (",nz,")."
      WRITE (*,*) "   Resetting same_res to 0."
      same_res = 0
    ENDIF

    WRITE(6,'(1x,2(a,f17.5))')                                          &
        'Original ctrlat1_old = ',ctrlat1_old,', ctrlon1_old = ',       &
        ctrlon1_old
    WRITE(6,'(1x,2(a,f17.5))')                                          &
        'Final ctrlat1        = ',ctrlat1    ,', ctrlon1     = ',ctrlon1
    WRITE(6,'(1x,2(a,f17.5))')                                          &
        'Final xctr1      = ',xctr1,  ', yctr1   = ',yctr1
    WRITE(6,'(a,i2/)')'same_res = ', same_res

    !wdt update: moved code block to above
!
!-----------------------------------------------------------------------
!
!  Perform spatial interpolations.
!
!-----------------------------------------------------------------------
!
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!    IF( nfile == 1) THEN
    IF( nfile == 1 .or. noutgrds > 0) THEN
       ! Need to do this for the first file only

!wdt update moved code block here
      xorig1 = xctr1 - (nx1-3)*dx1*0.5
      yorig1 = yctr1 - (ny1-3)*dy1*0.5
      zorig1 = zorig

      DO i=1,nx1
        x1(i) =xorig1+(i-2.0)*dx1
        xs1(i)=xorig1+(i-1.5)*dx1
      END DO

      DO j=1,ny1
        y1 (j) =yorig1+(j-2.0)*dy1
        ys1(j) =yorig1+(j-1.5)*dy1
      END DO

      ib= nint(x1(1)/dx)+1 
      jb= nint(y1(1)/dy)+1 

      DO i=1,nx-1
        xs(i)=0.5*(x(i)+x(i+1))
      END DO

      DO j=1,ny-1
        ys(j)=0.5*(y(j)+y(j+1))
      END DO

      DO k=1,nz1
        z1(k)=zorig1+(k-2.0)*dz1
      END DO

      IF( nfile == 1) THEN  ! Need to do this for the first file only
        xeps = 0.01*dx
        yeps = 0.01*dy

        PRINT*,'x (2),x (nx -1),y (2),y (ny -1) =',                       &
                x(2),x(nx-1),y(2),y(ny-1)
        PRINT*,'x1(2),x1(nx1-1),y1(2),y1(ny1-1) =',                       &
                x1(2),x1(nx1-1),y1(2),y1(ny1-1)

        IF(x1(    1) < x(   1)-xeps.OR.x1(nx1) > x(nx)+xeps.OR.           &
              y1(    1) < y(   1)-yeps.OR.y1(ny1) > y(ny)+yeps) THEN
          WRITE(6,'(3(/2x,a),/2x,2(a,f12.4),2(a,i5),2(a,f12.4),/2x,a)')   &
        'Sorry, at least part of your new grid is outside the border of', &
              'the original grid, please check input parameters',         &
              'dx1,dy1, nx1, ny1, xctr1 and yctr1. Currently,',           &
              'dx1=',dx1,', dy1=',dy1,', nx1=',nx1,', ny1=',ny1,          &
              ', xctr1=',xctr1,', yctr1=',yctr1,                          &
              'Job stopped in ARPSINTRP.'
          GOTO 600 !wdt
          STOP 1001
        END IF
      END IF

      IF( same_res == 0 ) THEN

        IF( intrphopt == 1) THEN ! Weights for 1st order interpolation

          DO i=1,nx1-1
            isx(i) = MAX(1, MIN(nx-2, INT((xs1(i)-xs(1))/dx)+1 ))
            wgtsx(i,1)= (xs(isx(i)+1)-xs1(i))/(xs(isx(i)+1)-xs(isx(i)))
          END DO

          DO j=1,ny1-1
            jsy(j) = MAX(1, MIN(ny-2, INT((ys1(j)-ys(1))/dy)+1 ))
            wgtsy(j,1)= (ys(jsy(j)+1)-ys1(j))/(ys(jsy(j)+1)-ys(jsy(j)))
          END DO

          DO i=1,nx1
            iux(i) = MAX(1, MIN(nx-1, INT((x1 (i)-x (1))/dx)+1 ))
            wgtux(i,1)= (x (iux(i)+1)-x1 (i))/(x (iux(i)+1)-x (iux(i)))
          END DO

          DO j=1,ny1
            jvy(j) = MAX(1, MIN(ny-1, INT((y1 (j)-y (1))/dy)+1 ))
            wgtvy(j,1)= (y (jvy(j)+1)-y1 (j))/(y (jvy(j)+1)-y (jvy(j)))
          END DO

        ELSE ! Weights for second order interpolation

          DO i=1,nx1-1
            isx(i) = MAX(2, MIN(nx-2, INT((xs1(i)-xs(1))/dx)+1 ))
            d=xs1(i)
            a=xs(isx(i)-1)
            b=xs(isx(i))
            c=xs(isx(i)+1)
            wgtsx(i,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtsx(i,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtsx(i,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO

          DO j=1,ny1-1
            jsy(j) = MAX(2, MIN(ny-2, INT((ys1(j)-ys(1))/dy)+1 ))
            d=ys1(j)
            a=ys(jsy(j)-1)
            b=ys(jsy(j))
            c=ys(jsy(j)+1)
            wgtsy(j,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtsy(j,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtsy(j,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO

          DO i=1,nx1
            iux(i) = MAX(2, MIN(nx-1, INT((x1 (i)-x (1))/dx)+1 ))
            d=x1(i)
            a=x(iux(i)-1)
            b=x(iux(i))
            c=x(iux(i)+1)
            wgtux(i,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtux(i,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtux(i,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO

          DO j=1,ny1
            jvy(j) = MAX(2, MIN(ny-1, INT((y1 (j)-y (1))/dy)+1 ))
            d=y1(j)
            a=y(jvy(j)-1)
            b=y(jvy(j))
            c=y(jvy(j)+1)
            wgtvy(j,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtvy(j,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtvy(j,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO

        END IF ! intrphopt = ?


        IF( ternopt1 == 0 ) THEN

          DO i=1,nx1-1
            DO j=1,ny1-1
              hterain1(i,j) = zrefsfc1
            END DO
          END DO

        ELSE IF( ternopt1 == 1 ) THEN  ! Bell-shaped mountain
!
!-----------------------------------------------------------------------
!
!  Define the bell-shaped mountain
!
!-----------------------------------------------------------------------
!
          pi2 = 1.5707963267949

          DO j=1,ny1-1
            DO i=1,nx1-1

              IF( mntwidy1 < 0.0 ) THEN   ! 2-d terrain in x-z plane.

                radnd = 1.0+((xs1(i)-mntctrx1)/mntwidx1)**2

              ELSE IF( mntwidx1 < 0.0 ) THEN ! 2-d terrain in y-z plane.

                radnd = 1.0+((ys1(j)-mntctry1)/mntwidy1)**2

              ELSE                             ! 3-d terrain

                radnd = 1.0+((xs1(i)-mntctrx1)/mntwidx1)**2             &
                       +((ys1(j)-mntctry1)/mntwidy1)**2
              END IF

              hterain1(i,j) = zrefsfc1 + hmount1/radnd

            END DO
          END DO

        ELSE IF( ternopt1 == 2 .or. ternopt1 == 4 ) THEN 
                                              ! Read from terrain data base

!
!-----------------------------------------------------------------------
!
!    Read in the terrain data.
!
!-----------------------------------------------------------------------
!
          lenstr = 132
          CALL strlnth( terndta1, lenstr)
          CALL readtrn( nx1,ny1,dx1,dy1, terndta1(1:lenstr),            &
                   mapproj,trulat1,trulat2,trulon,sclfct,ctrlat1,ctrlon1, &
                   hterain1 )

!
!-----------------------------------------------------------------------
!
!  Define the new grid terrain by interpolation
!
!-----------------------------------------------------------------------
!
!wdt update
        END IF

!wdt update
        IF( ternopt1 == 3 .or. ternopt1 == 4 ) THEN

!wdt update htrn1orig
          CALL intrpxy3d(hterain,nx,1,nx-1,ny,1,ny-1,1,1,1,             &
                         wgtsx,isx,wgtsy,jsy,intrphopt,                 &
                         htrn1orig,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

!wdt-williams 2002-01-17 GMB ntmerge
          !hterain1 = htrn1orig   ! uncomment if ternopt1=4 (ntmerge) dissabled
          IF (ternopt1 == 3) THEN
            hterain1 = htrn1orig
          ELSE 
            ntmergeinv = 1.d0/ntmerge
            DO j = 1,ny1
              DO i = 1,nx1
                idist = max(0,min(ntmerge,i-2,nx1-2-i,j-2,ny1-2-j)) 
                mfac = idist*ntmergeinv
                hterain1(i,j) = (1.d0-mfac)*htrn1orig(i,j)  &
                                + mfac*hterain1(i,j)
              END DO
            END DO
          END IF

        END IF

        npoint_below_ground = 0

        IF( ternopt1 /= 3) THEN
          CALL intrpxy3d(hterain,nx,1,nx-1,ny,1,ny-1,1,1,1,             &
                         wgtsx,isx,wgtsy,jsy,intrphopt,                 &
                         temx1y1z,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

          CALL a3dmax0(hterain,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,         &
                      amax,amin)
          WRITE(6,'(1x,2(a,e13.6))') 'htmin   = ', amin,', htmax    =',amax

          CALL a3dmax0(temx1y1z,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,    &
                      amax,amin)
          WRITE(6,'(1x,2(a,e13.6))') 'ht1min  = ', amin,', ht1max   =',amax

          CALL a3dmax0(hterain1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,    &
                      amax,amin)
          WRITE(6,'(1x,2(a,e13.6))') 'ht1min  = ', amin,', ht1max   =',amax

          npoint_below_ground = 0
          DO i=1,nx1-1
            DO j=1,ny1-1
              IF( hterain1(i,j) < temx1y1z(i,j,1)-1.0E-2) THEN
                npoint_below_ground = 1
                GO TO 2100
              END IF
            END DO
          END DO
          2100    CONTINUE
        END IF
!
!-----------------------------------------------------------------------
!
!  Set up a stretched vertical grid for the output grid.
!
!  For strhopt1=1, function y = a+b*x**3 is used to specify dz as a
!                              function of k.
!  For strhopt1=2, function y = c + a*tanh(b*x) is used to specify dz
!                              as a function of k.
!
!-----------------------------------------------------------------------
!

        IF ( strhopt1 == 0 ) THEN

          DO k=1,nz1
            zp1d1(k) = z1(k)
          END DO

        ELSE IF ( strhopt1 == 1 .OR.strhopt1 == 2 ) THEN

          za = zrefsfc1 + MAX(0.0, MIN(dlayer11, z1(nz1-2)-zrefsfc1 ))
          zb = za       + MAX(0.0, MIN(dlayer21, z1(nz1-1)-za      ))

          IF( dlayer11 >= (nz1-3)*dzmin1 ) THEN
            WRITE(6,'(/1x,a,f13.3,/a,f13.3,a,a)')                       &
                'Can not setup a vertical grid with uniform dz=',dzmin1, &
                ' over the depth of ',dlayer11,' please specify a smaller ', &
                'value of dlayer11. Program stopped in ARPSINTRP.'
            STOP 9105
          END IF

          CALL strhgrd(nz1,strhopt1,zrefsfc1,za,zb,z1(nz1-1),           &
                       dzmin1,strhtune1, zp1d1,dzp1d1)

        ELSE

          WRITE(6,'(1x,a,i3,a/)')                                       &
              'Invalid vertical grid stretching option, strhopt1 was ',strhopt1, &
              '. Program stopped in INIGRD.'
          STOP 9106

        END IF

!-----------------------------------------------------------------------
!
!  Physical height of computational grid defined as
!
!  Zp=(z-zrefsfc)*(Zm-hterain)/(Zm-zrefsfc)+hterain for z=<Zm.
!  ZP=z for z>Zm
!
!  where Zm the height at which the grid levels becomes flat.
!  Hm < Zm =< Ztop, hm is the height of mountain and Ztop the height
!  of model top.
!
!-----------------------------------------------------------------------
!
        DO k=nz1-1,2,-1
          IF(zp1d1(k) <= zflat1) THEN
            zflat11 = zp1d1(k)
            EXIT
          END IF
        END DO
        525   CONTINUE
        zflat11=MAX(MIN(z1(nz1-1),zflat11),zrefsfc1)

        DO k=2,nz1-1

          IF(zp1d1(k) > zflat11) THEN
            DO j=1,ny1-1
              DO i=1,nx1-1
                zp1(i,j,k)=zp1d1(k)
              END DO
            END DO
          ELSE
            DO j=1,ny1-1
              DO i=1,nx1-1
                zp1(i,j,k)=(zp1d1(k)-zrefsfc1)*(zflat11-hterain1(i,j))  &
                           /(zflat11-zrefsfc1)+hterain1(i,j)
              END DO
            END DO
          END IF

        END DO

        DO j=1,ny1-1
          DO i=1,nx1-1
            zp1(i,j,2)=hterain1(i,j)
            zp1(i,j,1)=2.0*zp1(i,j,2)-zp1(i,j,3)
            zp1(i,j,nz1)=2.0*zp1(i,j,nz1-1)-zp1(i,j,nz1-2)
          END DO
        END DO
!
        CALL jacob(nx1,ny1,nz1,x1,y1,z1,zp1,j11,j21,j31)

      END IF ! same_res = 0?

    END IF  ! IF (nfile.eq.1)
!
!-----------------------------------------------------------------------
!
!  Calculate total fields from that for base state and perturbations
!
!-----------------------------------------------------------------------
!
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          u (i,j,k)=uprt (i,j,k)+ubar (i,j,k)
          v (i,j,k)=vprt (i,j,k)+vbar (i,j,k)
          qv(i,j,k)=qvprt(i,j,k)+qvbar(i,j,k)
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Start interpolation or sampling in the case of same_res=1.
!
!-----------------------------------------------------------------------
!
    IF( same_res == 1 ) THEN

      DO k=1,nz1
        DO j=1,ny1
          DO i=1,nx1
            u1    (i,j,k)=u    (i+ib,j+jb,k)
            v1    (i,j,k)=v    (i+ib,j+jb,k)
            w1    (i,j,k)=w    (i+ib,j+jb,k)
            ptprt1(i,j,k)=ptprt(i+ib,j+jb,k)
            pprt1 (i,j,k)=pprt (i+ib,j+jb,k)
            qv1   (i,j,k)=qv   (i+ib,j+jb,k)
            qc1   (i,j,k)=qc   (i+ib,j+jb,k)
            qr1   (i,j,k)=qr   (i+ib,j+jb,k)
            qi1   (i,j,k)=qi   (i+ib,j+jb,k)
            qs1   (i,j,k)=qs   (i+ib,j+jb,k)
            qh1   (i,j,k)=qh   (i+ib,j+jb,k)
            tke1  (i,j,k)=tke  (i+ib,j+jb,k)
            kmh1  (i,j,k)=kmh  (i+ib,j+jb,k)
            kmv1  (i,j,k)=kmv  (i+ib,j+jb,k)
            radfrc1(i,j,k)=radfrc(i+ib,j+jb,k)
            rhobar1(i,j,k)=rhobar(i+ib,j+jb,k)
          END DO
        END DO
      END DO

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!      IF( nfile == 1) THEN
      IF( nfile == 1 .or. noutgrds > 0) THEN

        DO k=1,nz1
          DO i=1,nx1
            DO j=1,ny1
              ubar1 (i,j,k)=ubar (i+ib,j+jb,k)
              vbar1 (i,j,k)=vbar (i+ib,j+jb,k)
              ptbar1(i,j,k)=ptbar(i+ib,j+jb,k)
              pbar1 (i,j,k)=pbar (i+ib,j+jb,k)
              qvbar1(i,j,k)=qvbar(i+ib,j+jb,k)
              zp1   (i,j,k)=zp   (i+ib,j+jb,k)
            END DO
          END DO
        END DO

        CALL jacob(nx1,ny1,nz1,x1,y1,z1,zp1,j11,j21,j31)

      END IF

      DO i=1,nx1
        DO j=1,ny1
          hterain1(i,j)=hterain(i+ib,j+jb)

          soiltyp1 (i,j,:)=soiltyp (i+ib,j+jb,:)
          stypfrct1(i,j,:)=stypfrct(i+ib,j+jb,:)

          vegtyp1(i,j)=vegtyp(i+ib,j+jb)
          lai1   (i,j)=lai   (i+ib,j+jb)
          roufns1(i,j)=roufns(i+ib,j+jb)
          veg1   (i,j)=veg   (i+ib,j+jb)

          tsfc1    (i,j,:)=tsfc    (i+ib,j+jb,:)
          tsoil1   (i,j,:)=tsoil   (i+ib,j+jb,:)
          wetsfc1  (i,j,:)=wetsfc  (i+ib,j+jb,:)
          wetdp1   (i,j,:)=wetdp   (i+ib,j+jb,:)
          wetcanp1 (i,j,:)=wetcanp (i+ib,j+jb,:)
          snowdpth1(i,j)=snowdpth(i+ib,j+jb)

          raing1  (i,j)=raing(i+ib,j+jb)
          rainc1  (i,j)=rainc(i+ib,j+jb)

          prcrate1(i,j,1:4)=prcrate(i+ib,j+jb,1:4)

          radsw1 (i,j)=radsw (i+ib,j+jb)
          rnflx1 (i,j)=rnflx (i+ib,j+jb)
          usflx1 (i,j)=usflx (i+ib,j+jb)
          vsflx1 (i,j)=vsflx (i+ib,j+jb)
          ptsflx1(i,j)=ptsflx(i+ib,j+jb)
          qvsflx1(i,j)=qvsflx(i+ib,j+jb)
        END DO
      END DO

      GO TO 800

    END IF
!
!-----------------------------------------------------------------------
!
!  Intepolate the 2-D arrays to the new grid
!
!-----------------------------------------------------------------------
!

    CALL intrpxy3d(rainc,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   rainc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(raing,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   raing1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(prcrate,nx,1,nx-1,ny,1,ny-1,4,1,4,                   &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   prcrate1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(radsw,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   radsw1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(rnflx,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   rnflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(usflx,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   usflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(vsflx,nx,1,nx-1,ny,1,ny-1,1,1,1,                     &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   vsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(ptsflx,nx,1,nx-1,ny,1,ny-1,1,1,1,                    &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   ptsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(qvsflx,nx,1,nx-1,ny,1,ny-1,1,1,1,                    &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   qvsflx1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(lai,nx,1,nx-1,ny,1,ny-1,1,1,1,                       &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   lai1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(roufns,nx,1,nx-1,ny,1,ny-1,1,1,1,                    &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   roufns1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(veg,nx,1,nx-1,ny,1,ny-1,1,1,1,                       &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   veg1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
    DO j = 1,ny1
      DO i = 1,nx1
        ix(i,j) = isx(i)
        jy(i,j) = jsy(j)
        xw(i,j) =  &
           MAX(0.,MIN(1.,(xs(isx(i)+1)-xs1(i))/(xs(isx(i)+1)-xs(isx(i)))))
        yw(i,j) =  &
           MAX(0.,MIN(1.,(ys(jsy(j)+1)-ys1(j))/(ys(jsy(j)+1)-ys(jsy(j)))))
      END DO
    END DO
    CALL intrp_soil(nx,ny,nx1,ny1,nstypin,nstyp1,xw,yw,ix,jy, &
                    tsfc,tsoil,wetsfc,wetdp,wetcanp,  &
                    soiltyp,stypfrct,vegtyp, &
                    tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,  &
                    soiltyp1,stypfrct1,vegtyp1)
!    CALL intrpxy3d(tsfc,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1,          &
!                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
!                   tsfc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
!    CALL intrpxy3d(tsoil,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1,         &
!                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
!                   tsoil1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
!    CALL intrpxy3d(wetsfc,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1,        &
!                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
!                   wetsfc1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
!    CALL intrpxy3d(wetdp,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1,         &
!                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
!                   wetdp1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)
!
!    CALL intrpxy3d(wetcanp,nx,1,nx-1,ny,1,ny-1,nstyp+1,1,nstyp+1,       &
!                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
!                   wetcanp1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

    CALL intrpxy3d(snowdpth,nx,1,nx-1,ny,1,ny-1,1,1,1,                  &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   snowdpth1,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!    DO i=1,nx1-1
!      DO j=1,ny1-1
!        vegtyp1 (i,j) = vegtyp (isx(i),jsy(j))
!      END DO
!    END DO
!
!    DO iss=1,nstyp
!      DO j=1,ny1-1
!        DO i=1,nx1-1
!          soiltyp1 (i,j,iss) = soiltyp (isx(i),jsy(j),iss)
!          stypfrct1(i,j,iss) = stypfrct(isx(i),jsy(j),iss)
!        END DO
!      END DO
!    END DO

!
!-----------------------------------------------------------------------
!
!  Intepolate 3D arrays
!
!-----------------------------------------------------------------------
!
    IF( z_or_p == 2) THEN  ! This situation has not been tested since
                           ! last major change to this program
      DO k=1,nz
        DO j=1,ny-1
          DO i=1,nx-1
            zp(i,j,k) = -ALOG( (pbar(i,j,k)+pprt(i,j,k))*1.0E-5 )
          END DO
        END DO
      END DO

      DO k=1,nz1
        DO j=1,ny1-1
          DO i=1,nx1-1
            zp1(i,j,k) = -ALOG( plevels(k)*1.0E-3 )
          END DO
        END DO
      END DO
    END IF


    IF( npoint_below_ground /= 0.AND.z_or_p == 1.AND.                   &
          (bglopt == 4 .OR. bglopt == 5)) THEN
          ! Reconstruct horizontally homogeneous base state
      redo_base_state = 1
    ELSE
      redo_base_state = 0
    END IF

    IF( redo_base_state /= 0 ) THEN

      zpmin=0.5*(zp(2,2,    2)+zp(2,2,3))
      zpmax=0.5*(zp(2,2,nz)+zp(2,2,nz-1))
      DO i=2,nx-2
        DO j=2,ny-2
          zpmin=MIN(zpmin,0.5*(zp(i,j, 2)+zp(i,j,   3)))
          zpmax=MAX(zpmax,0.5*(zp(i,j,nz)+zp(i,j,nz-1)))
        END DO
      END DO

      zpmin1=MIN(zpmin1,0.5*(zp1(1,1,1)+zp1(1,1,2)))
      DO i=2,nx1-2
        DO j=2,ny1-2
          zpmin1=MIN(zpmin1,0.5*(zp1(i,j,1)+zp1(i,j,2)))
        END DO
      END DO

      snddelz = (zpmax-zpmin)/(lvlprof-2)
      DO k=1,lvlprof
        zsnd(k)= zpmin+(k-2)*snddelz
      END DO

      IF (bglopt == 4) THEN

        DO k=2,lvlprof
          temz1d3(k)=0.0
          temz1d4(k)=0.0
          temz1d5(k)=0.0
          temz1d6(k)=0.0
          temz1d7(k)=0.0
        END DO
        pbartop = 0.0

        DO i=2,nx-2
          DO j=2,ny-2

            DO k=2,nz-1
              temz1d1(k)=ptbar(i,j,k)
              temz1d2(k)=0.5*(zp(i,j,k)+zp(i,j,k+1))
            END DO
            CALL inte1d(temz1d1(2),temz1d2(2),nz-2,                     &
                      ptsnd(2),zsnd(2),lvlprof-1)

            DO k=2,nz-1
              temz1d1(k)=qvbar(i,j,k)
            END DO
            CALL inte1d(temz1d1(2),temz1d2(2),nz-2,                     &
                      qvsnd(2),zsnd(2),lvlprof-1)

            DO k=2,nz-1
              temz1d1(k)=(ubar(i,j,k)+ubar(i+1,j,k))*0.5
            END DO
            CALL inte1d(temz1d1(2),temz1d2(2),nz-2,                     &
                      ubsnd(2),zsnd(2),lvlprof-1)

            DO k=2,nz-1
              temz1d1(k)=(vbar(i,j,k)+vbar(i,j+1,k))*0.5
            END DO
            CALL inte1d(temz1d1(2),temz1d2(2),nz-2,                     &
                      vbsnd(2),zsnd(2),lvlprof-1)

            DO k=2,nz-1
              temz1d1(k)=ALOG(pbar(i,j,k))
            END DO
            CALL inte1d(temz1d1(2),temz1d2(2),nz-2,                     &
                      tem,zsnd(lvlprof-1),1)

            pbartop = pbartop+tem

            DO k=2,lvlprof
              IF(zsnd(k) >= temz1d2(2)-1.0E-5) THEN
                 ! Add only points above ground
                temz1d3(k)=temz1d3(k)+1
                temz1d4(k)=temz1d4(k)+ptsnd(k)
                temz1d5(k)=temz1d5(k)+qvsnd(k)
                temz1d6(k)=temz1d6(k)+ubsnd(k)
                temz1d7(k)=temz1d7(k)+vbsnd(k)
              END IF
            END DO

          END DO
        END DO
!
!  Obtain averages
!
        DO k=2,lvlprof
          ptsnd(k)=temz1d4(k)/temz1d3(k)
          qvsnd(k)=temz1d5(k)/temz1d3(k)
          ubsnd(k)=temz1d6(k)/temz1d3(k)
          vbsnd(k)=temz1d7(k)/temz1d3(k)
        END DO

        pbartop=EXP(pbartop/temz1d3(lvlprof-1))
!
!  Integrate hydrostatic equation down from top
!
        DO k=2,lvlprof
          temz1d1(k) = ptsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
        END DO

        psnd(lvlprof-1) = ( pbartop/p0 )**(rddcp)

        DO k=lvlprof-2,2,-1
          psnd(k)=psnd(k+1)-g*(zsnd(k)-zsnd(k+1))                       &
                       /(cp*0.5*(temz1d1(k)+temz1d1(k+1)))
        END DO
        DO k=lvlprof,lvlprof
          psnd(k)=psnd(k-1)-g*(zsnd(k)-zsnd(k-1))                       &
                       /(cp*0.5*(temz1d1(k)+temz1d1(k-1)))
        END DO

        DO k = 2,lvlprof
          psnd(k) = psnd(k)**(1/rddcp) * p0
        END DO

        IF( zpmin1 < zpmin ) THEN  ! Extend 1d sounding to zpmin1 using
                                   ! constant temperature lapse rate

          dtdz = -0.0068  ! Minus lapse rate
          zsnd(1)=zpmin1
          delz = zpmin1 - zpmin

          tbark1 = ptsnd(2) * (psnd(2)/p0)**rddcp
          tbark  = tbark1+dtdz*delz

          lnpbar = ALOG(psnd(2))-g*2.0*delz/(rd*(tbark1+tbark))
          psnd(1) = EXP(lnpbar)

          ptsnd(1)=tbark*(p0/psnd(1))**rddcp
          qvsnd(1)=qvsnd(2)
          ubsnd(1)=ubsnd(2)
          vbsnd(1)=vbsnd(2)

        ELSE  ! Set for safty, won't be used.

          psnd(1)=psnd(2)
          ptsnd(1)=ptsnd(2)
          qvsnd(1)=qvsnd(2)
          ubsnd(1)=ubsnd(2)
          vbsnd(1)=vbsnd(2)

        END IF

      ELSE ! bglopt=5
!wdt update
!
!  Use extmnsnd (from ext2arps) to get average sounding
!
        ! tem1 - zp at zone center
        DO k=1,nz
          DO j=1,ny
            DO i=1,nx
              tem1(i,j,k) = 0.5*(zp(i,j,k)+zp(i,j,k+1))
            END DO
          END DO
        END DO

        ! tem2 = log(pbar)
        ! tem3 = t
        ! tem4 = rhstar
        DO k=1,nz
          DO j=1,ny
            DO i=1,nx
              tem2(i,j,k) = ALOG(pbar(i,j,k))                      ! log pressure
              tem3(i,j,k) = (ptprt(i,j,k)+ptbar(i,j,k))                 &
                         *(((pprt(i,j,k)+pbar(i,j,k))/p0)**rddcp)  ! temperature
              qvsat=f_qvsat( pprt(i,j,k)+pbar(i,j,k), tem3(i,j,k) )
              tem4(i,j,k)=SQRT(AMAX1(0.,(rhmax-(qv(i,j,k)/qvsat))))
            END DO
          END DO
        END DO

        DO j=1,ny-1
          DO i=1,nx-1
            xs_2d(i,j) = xs(i)
            ys_2d(i,j) = ys(j)
          END DO
        END DO

        i1mn = MAX(1,INT(xorig1/dx1-1.))
        i1mx = MIN(nx-1,i1mn+nx1+1)
        j1mn = MAX(1,INT(yorig1/dy1-1.))
        j1mx = MIN(ny-1,j1mn+ny1+1)

!wdt update
        CALL extmnsnd(nx1,ny1,nx,ny,nz,lvlprof,                         &
                        i1mn,i1mx,j1mn,j1mx,2,nz-1,                     &
                        xs1,ys1,xs_2d,ys_2d,                            &
                        tem1,tem2,tem3,tem4,                            &
                        u,v,   & ! note: not vital that u & v be staggered
                        zsnd,temz1d1,psnd,temz1d2,ptsnd,temz1d3,        &
                        qvsnd,temz1d4,ubsnd,vbsnd)

      END IF

      WRITE(6,'(/a/)')'    Reconstructed horizontal average sounding:'
      DO k = lvlprof,1,-1
        WRITE(6,'(a,i4,4f10.2,f10.6)') 'k,z,p,pt,t,qv=',                &
                k,zsnd(k),psnd(k),ptsnd(k),                             &
                ptsnd(k)*(psnd(k)/p0)**rddcp,qvsnd(k)
      END DO

    END IF
!
!-----------------------------------------------------------------------
!
!  Intepolate 3D arrays at u-point
!
!-----------------------------------------------------------------------
!
    CALL avgz(zp1, 0 , nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1, tem11 )
    CALL avgsu(tem11,nx1,ny1,nz1,1,ny1-1,1,nz1-1, tem21, tem11)

    CALL avgz(zp, 0 , nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem1 )
    CALL avgsu(tem1,nx,ny,nz,1,ny-1,1,nz-1, tem2, tem1)

    CALL intrpxy3d(tem2,nx,1,nx,ny,1,ny-1,nz,1,nz-1,                    &
                   wgtux,iux,wgtsy,jsy,intrphopt,                       &
                   temx1y1zb,nx1,1,nx1, ny1,1,ny1-1, temx1yz)

    DO i=1,nx1
      DO j=1,ny1-1
        k0 = nz-2
        DO k=nz1-1,1,-1
          z0 = tem21(i,j,k)
          DO kk=k0,1,-1
            IF(z0 >= temx1y1zb(i,j,kk)) THEN
              k0 = kk
              EXIT
            END IF
          END DO
          1115      CONTINUE
          IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
          IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
        END DO
      END DO
    END DO

    IF( intrpvopt == 1) THEN
      DO i=1,nx1
        DO j=1,ny1-1
          DO k=1,nz1-1
            wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/    &
                     (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
          END DO
        END DO
      END DO
    ELSE
      DO i=1,nx1
        DO j=1,ny1-1
          DO k=1,nz1-1
            d=tem21(i,j,k)
            a=temx1y1zb(i,j,kz(i,j,k)-1)
            b=temx1y1zb(i,j,kz(i,j,k))
            c=temx1y1zb(i,j,kz(i,j,k)+1)
            wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO
        END DO
      END DO
    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!    IF( nfile == 1) THEN
    IF(nfile == 1 .or. noutgrds > 0) THEN

      IF( redo_base_state == 0 ) THEN ! get base-state from original arrays

        CALL intrpxyz3d(ubar,nx,1,nx,ny,1,ny-1,nz,1,nz-1,               &
                        wgtux,iux,wgtsy,jsy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        ubar1,nx1,1,nx1, ny1,1,ny1-1,nz1,1,nz1-1,       &
                        temx1yz,temx1y1z)
        DO i=1,nx1
          DO j=1,ny1-1
            DO k=1,nz1-1
              IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
                IF(bglopt == 2) THEN
                  ubar1(i,j,k)=temx1y1z(i,j,2)
                ELSE IF(bglopt == 3) THEN
                  ubar1(i,j,k)=0.0
                END IF
              END IF
            END DO
          END DO
        END DO

      ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings

        DO k=1,nz1-1
          DO j=1,ny1-1
            DO i=1,nx1
              kk = MAX(1,MIN(lvlprof-1,                                 &
                   INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
              IF(tem21(i,j,k) < zsnd(2)) kk = 1
              tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
              ubar1(i,j,k)= ubsnd(kk)+(ubsnd(kk+1)-ubsnd(kk))*tem
            END DO
          END DO
        END DO

      END IF

    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
                           ! same as input surface wind (by matching tem21
                           ! to the input height, temx1y1zb, near the surface).
      DO i=1,nx1
        DO j=1,ny1-1
          hght0 = tem21(i,j,2)
          dhght = hght0 - temx1y1zb(i,j,2)
          IF (dhght > 0) THEN
            DO k=1,nz1-1
              fac = MAX(0., tem21(i,j,k) - hght0)/aghght
              IF (fac > 1) EXIT
              tem21(i,j,k) = fac*fac*tem21(i,j,k)  &
                             + (1.-fac*fac)*(tem21(i,j,k)-dhght)
            END DO
          ENDIF
        END DO
      END DO

      DO i=1,nx1
        DO j=1,ny1-1
          k0 = nz-2
          DO k=nz1-1,1,-1
            z0 = tem21(i,j,k)
            DO kk=k0,1,-1
              IF(z0 >= temx1y1zb(i,j,kk)) THEN
                k0 = kk
                EXIT
              END IF
            END DO
            IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
            IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
          END DO
        END DO
      END DO

      IF( intrpvopt == 1) THEN
        DO i=1,nx1
          DO j=1,ny1-1
            DO k=1,nz1-1
              wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/    &
                       (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
            END DO
          END DO
        END DO
      ELSE
        DO i=1,nx1
          DO j=1,ny1-1
            DO k=1,nz1-1
              d=tem21(i,j,k)
              a=temx1y1zb(i,j,kz(i,j,k)-1)
              b=temx1y1zb(i,j,kz(i,j,k))
              c=temx1y1zb(i,j,kz(i,j,k)+1)
              wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
              wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
              wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
            END DO
          END DO
        END DO
      END IF

    END IF

    CALL intrpxyz3d(u,nx,1,nx,ny,1,ny-1,nz,1,nz-1,                      &
                    wgtux,iux,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    u1,nx1,1,nx1, ny1,1,ny1-1,nz1,1,nz1-1,              &
                    temx1yz,temx1y1z)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground

      DO i=1,nx1
        DO j=1,ny1-1
          DO k=nz1-2,1,-1
            IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
              IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
                u1   (i,j,k)=temx1y1z(i,j,2)
              ELSE IF(bglopt == 3) THEN
                u1   (i,j,k)=misvalue
              END IF
            END IF
          END DO
        END DO
      END DO

    END IF
!
!-----------------------------------------------------------------------
!
!  Intepolate 3D arrays at v-point
!
!-----------------------------------------------------------------------
!
    CALL avgz(zp1, 0, nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1, tem11 )
    CALL avgsv(tem11,nx1,ny1,nz1,1,nx1-1,1,nz1-1, tem21, tem11)

    CALL avgz(zp, 0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, tem1 )
    CALL avgsv(tem1,nx,ny,nz,1,nx-1,1,nz-1, tem2, tem1)

    CALL intrpxy3d(tem2,nx,1,nx-1,ny,1,ny,nz,1,nz-1,                    &
                   wgtsx,isx,wgtvy,jvy,intrphopt,                       &
                   temx1y1zb,nx1,1,nx1-1, ny1,1,ny1, temx1yz)

    DO i=1,nx1-1
      DO j=1,ny1

        k0=nz-2
        DO k=nz1-1,1,-1
          z0 = tem21(i,j,k)
          DO kk=k0,1,-1
            IF(z0 >= temx1y1zb(i,j,kk)) THEN
              k0 = kk
              EXIT
            END IF
          END DO
          1215      CONTINUE
          IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
          IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
        END DO

      END DO
    END DO

    IF( intrpvopt == 1) THEN
      DO i=1,nx1-1
        DO j=1,ny1
          DO k=1,nz1-1
            wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/    &
                      (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
          END DO
        END DO
      END DO
    ELSE
      DO i=1,nx1-1
        DO j=1,ny1
          DO k=1,nz1-1
            d=tem21(i,j,k)
            a=temx1y1zb(i,j,kz(i,j,k)-1)
            b=temx1y1zb(i,j,kz(i,j,k))
            c=temx1y1zb(i,j,kz(i,j,k)+1)
            wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO
        END DO
      END DO
    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!    IF( nfile == 1) THEN
    IF(nfile == 1 .or. noutgrds > 0) THEN

      IF( redo_base_state == 0 ) THEN ! get base-state from original arrays

        CALL intrpxyz3d(vbar,nx,1,nx-1,ny,1,ny,nz,1,nz-1,               &
                        wgtsx,isx,wgtvy,jvy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        vbar1,nx1,1,nx1-1, ny1,1,ny1,nz1,1,nz1-1,       &
                        temx1yz,temx1y1z)
        DO i=1,nx1-1
          DO j=1,ny1
            DO k=nz1-2,1,-1
              IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
                IF(bglopt == 2) THEN
                  vbar1(i,j,k)=temx1y1z(i,j,2)
                ELSE IF(bglopt == 3) THEN
                  vbar1(i,j,k)=0.0
                END IF
              END IF
            END DO
          END DO
        END DO

      ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings

        DO k=1,nz1-1
          DO i=1,nx1-1
            DO j=1,ny1
              kk = MAX(1,MIN(lvlprof-1,                                 &
                   INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
              IF(tem21(i,j,k) < zsnd(2)) kk = 1
              tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
              vbar1(i,j,k)= vbsnd(kk)+(vbsnd(kk+1)-vbsnd(kk))*tem
            END DO
          END DO
        END DO

      END IF

    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
                           ! same as input surface wind (by matching tem21
                           ! to the input height, temx1y1zb, near the surface).
      DO i=1,nx1-1
        DO j=1,ny1
          hght0 = tem21(i,j,2)
          dhght = hght0 - temx1y1zb(i,j,2)
          IF (dhght > 0) THEN
            DO k=1,nz1-1
              fac = MAX(0., tem21(i,j,k) - hght0)/aghght
              IF (fac > 1) EXIT
              tem21(i,j,k) = fac*fac*tem21(i,j,k)  &
                             + (1.-fac*fac)*(tem21(i,j,k)-dhght)
            END DO
          ENDIF
        END DO
      END DO

      DO i=1,nx1-1
        DO j=1,ny1

          k0=nz-2
          DO k=nz1-1,1,-1
            z0 = tem21(i,j,k)
            DO kk=k0,1,-1
              IF(z0 >= temx1y1zb(i,j,kk)) THEN
                k0 = kk
                EXIT
              END IF
            END DO
            IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
            IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
          END DO

        END DO
      END DO

      IF( intrpvopt == 1) THEN
        DO i=1,nx1-1
          DO j=1,ny1
            DO k=1,nz1-1
              wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/    &
                        (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
            END DO
          END DO
        END DO
      ELSE
        DO i=1,nx1-1
          DO j=1,ny1
            DO k=1,nz1-1
              d=tem21(i,j,k)
              a=temx1y1zb(i,j,kz(i,j,k)-1)
              b=temx1y1zb(i,j,kz(i,j,k))
              c=temx1y1zb(i,j,kz(i,j,k)+1)
              wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
              wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
              wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
            END DO
          END DO
        END DO
      END IF

    END IF

    CALL intrpxyz3d(v,nx,1,nx-1,ny,1,ny,nz,1,nz-1,                      &
                    wgtsx,isx,wgtvy,jvy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    v1,nx1,1,nx1-1, ny1,1,ny1,nz1,1,nz1-1,              &
                    temx1yz,temx1y1z)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground

      DO i=1,nx1-1
        DO j=1,ny1
          DO k=nz1-2,1,-1
            IF(tem21(i,j,k) < temx1y1zb(i,j,2)) THEN ! below input grid ground
              IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
                v1   (i,j,k)=temx1y1z(i,j,2)
              ELSE IF(bglopt == 3) THEN
                v1   (i,j,k)=misvalue
              END IF
            END IF
          END DO
        END DO
      END DO

    END IF
!
!-----------------------------------------------------------------------
!
!  Intepolate 3D arrays at w-point
!
!-----------------------------------------------------------------------
!
    CALL intrpxy3d(zp,nx,1,nx-1,ny,1,ny-1,nz,1,nz,                      &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   temx1y1zb,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)

!wdt update: use tem21 instead of zp1 below
    tem21 = zp1

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF (ntagopt == 1) THEN ! adjust tem21 (heights) so that surface wind is
                          ! same as input surface wind (by matching tem21
                          ! to the input height, temx1y1zb, near the surface).
      DO i=1,nx1-1
        DO j=1,ny1-1
          hght0 = tem21(i,j,2)
          dhght = hght0 - temx1y1zb(i,j,2)
          IF (dhght > 0) THEN
            DO k=1,nz1-1
              fac = MAX(0., tem21(i,j,k) - hght0)/aghght
              IF (fac > 1) EXIT
              tem21(i,j,k) = fac*fac*tem21(i,j,k)  &
                             + (1.-fac*fac)*(tem21(i,j,k)-dhght)
            END DO
          ENDIF
        END DO
      END DO

    END IF

    DO i=1,nx1-1
      DO j=1,ny1-1
        k0 = nz-1
        DO k=nz1,1,-1
          z0 = tem21(i,j,k)
          DO kk=k0,1,-1
            IF(z0 >= temx1y1zb(i,j,kk)) THEN
              k0 = kk
              EXIT
            END IF
          END DO
          1315      CONTINUE

          IF(intrpvopt == 1) kz(i,j,k)=MAX(2,MIN(k0,nz-1))
          IF(intrpvopt == 2) kz(i,j,k)=MAX(3,MIN(k0,nz-1))
        END DO

      END DO
    END DO

    IF( intrpvopt == 1) THEN
      DO i=1,nx1-1
        DO j=1,ny1-1
          DO k=1,nz1
            wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/      &
                      (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
          END DO
        END DO
      END DO
    ELSE
      DO i=1,nx1-1
        DO j=1,ny1-1
          DO k=1,nz1
            d=tem21(i,j,k)
            a=temx1y1zb(i,j,kz(i,j,k)-1)
            b=temx1y1zb(i,j,kz(i,j,k))
            c=temx1y1zb(i,j,kz(i,j,k)+1)
            wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO
        END DO
      END DO
    END IF

    CALL intrpxyz3d(w,nx,1,nx-1,ny,1,ny-1,nz,1,nz,                      &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    w1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1,              &
                    temx1yz,temx1y1z)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.: ntagopt,aghght
    IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground

      DO i=1,nx1-1
        DO j=1,ny1-1
          DO k=nz1-1,1,-1
            IF(zp1(i,j,k) < temx1y1zb(i,j,2)) THEN ! Blow input grid ground
              IF(bglopt == 2 .OR. bglopt == 4 .OR. bglopt == 5) THEN
                w1(i,j,k)=temx1y1z(i,j,2)
              ELSE IF(bglopt == 3) THEN
                w1(i,j,k)=misvalue
              END IF
            END IF
          END DO
        END DO
      END DO

    END IF

!-----------------------------------------------------------------------
!
!  Interpolate 3D arrays at scalar point
!
!-----------------------------------------------------------------------

    CALL avgz(zp1,0,nx1,ny1,nz1,1,nx1-1,1,ny1-1,1,nz1-1,tem21)

    CALL avgz(zp,0,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem2)

    CALL intrpxy3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                  &
                   wgtsx,isx,wgtsy,jsy,intrphopt,                       &
                   temx1y1zb,nx1,1,nx1-1, ny1,1,ny1-1, temx1yz)


    DO i=1,nx1-1
      DO j=1,ny1-1
        htrnx1y1(i,j)=temx1y1zb(i,j,2) ! Old terrain on new grid at scalar level
                                       ! e.g., z of first grid level AGL.
        DO k=1,nz1-1
          IF(tem21(i,j,k) >= htrnx1y1(i,j)) GO TO 1350
        END DO
        1350    CONTINUE
        ktrnx1y1(i,j)= k-1 ! The first scale-level below the first scale level
                           ! above input grid ground
      END DO
    END DO

    DO i=1,nx1-1
      DO j=1,ny1-1

        k0=nz-2
        DO k=nz1-1,1,-1
          z0 = tem21(i,j,k)
          DO kk=k0,1,-1
            IF(z0 >= temx1y1zb(i,j,kk) ) THEN
              k0 = kk
              EXIT
            END IF
          END DO
          1415      CONTINUE
          IF(intrpvopt == 1) kz(i,j,k)=MAX(1,MIN(k0,nz-2))
          IF(intrpvopt == 2) kz(i,j,k)=MAX(2,MIN(k0,nz-2))
        END DO

      END DO
    END DO

    IF( intrpvopt == 1) THEN
      DO i=1,nx1-1
        DO j=1,ny1-1
          DO k=1,nz1-1
            wgtz(i,j,k,1)=(temx1y1zb(i,j,kz(i,j,k)+1)-tem21(i,j,k))/    &
                      (temx1y1zb(i,j,kz(i,j,k)+1)-temx1y1zb(i,j,kz(i,j,k)))
          END DO
        END DO
      END DO
    ELSE
      DO i=1,nx1-1
        DO j=1,ny1-1
          DO k=1,nz1-1
            d=tem21(i,j,k)
            a=temx1y1zb(i,j,kz(i,j,k)-1)
            b=temx1y1zb(i,j,kz(i,j,k))
            c=temx1y1zb(i,j,kz(i,j,k)+1)
            wgtz(i,j,k,1)= (d-b)*(d-c)/((a-b)*(a-c))
            wgtz(i,j,k,2)= (d-a)*(d-c)/((b-a)*(b-c))
            wgtz(i,j,k,3)= (d-a)*(d-b)/((c-a)*(c-b))
          END DO
        END DO
      END DO
    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!    IF( nfile == 1) THEN
    IF(nfile == 1 .or. noutgrds > 0) THEN

      IF( redo_base_state == 0) THEN ! Get base-state from original arrays

        DO k=1,nz-1
          DO j=1,ny-1
            DO i=1,nx-1
              tem1(i,j,k)=ALOG(pbar(i,j,k))
            END DO
          END DO
        END DO

        CALL intrpxyz3d(tem1,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,             &
                        wgtsx,isx,wgtsy,jsy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        pbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,     &
                        temx1yz,temx1y1z)

        DO k=1,nz1-1
          DO j=1,ny1-1
            DO i=1,nx1-1
              pbar1(i,j,k)=EXP(pbar1(i,j,k))
            END DO
          END DO
        END DO

        DO i=1,nx1-1
          DO j=1,ny1-1
            pbsfc(i,j)= EXP(temx1y1z(i,j,2))
          END DO
        END DO

        CALL intrpxyz3d(ptbar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,            &
                        wgtsx,isx,wgtsy,jsy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        ptbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,    &
                        temx1yz,temx1y1z)

        DO i=1,nx1-1
          DO j=1,ny1-1
            ptbsfc(i,j)= temx1y1z(i,j,2)
          END DO
        END DO


        CALL intrpxyz3d(qvbar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,            &
                        wgtsx,isx,wgtsy,jsy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        qvbar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,    &
                        temx1yz,temx1y1z)

        DO i=1,nx1-1
          DO j=1,ny1-1
            qvbsfc(i,j)= temx1y1z(i,j,2)
          END DO
        END DO

        CALL intrpxyz3d(rhobar,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,           &
                        wgtsx,isx,wgtsy,jsy,wgtz,kz,                    &
                        intrphopt,intrpvopt,                            &
                        rhobar1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,   &
                        temx1yz,temx1y1z)

      ELSE IF( redo_base_state == 1 ) THEN ! reconstruct base-state from 1d soundings

        DO k=1,nz1-1
          DO i=1,nx1-1
            DO j=1,ny1-1

              kk = MAX(1,MIN(lvlprof-1,                                 &
                   INT((tem21(i,j,k)-zsnd(2))/snddelz)+2))
              IF(tem21(i,j,k) < zsnd(2)) kk = 1

              tem = (tem21(i,j,k)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
              ptbar1(i,j,k)= ptsnd(kk)+(ptsnd(kk+1)-ptsnd(kk))*tem
              qvbar1(i,j,k)= qvsnd(kk)+(qvsnd(kk+1)-qvsnd(kk))*tem
              tema = ALOG(psnd (kk) )
              pbar1 (i,j,k)=EXP(tema+(ALOG(psnd(kk+1))-tema)*tem)
              rhobar1(i,j,k)=pbar1(i,j,k)/                              &
                            (rd*ptbar1(i,j,k)*(pbar1(i,j,k)/p0)**rddcp)

            END DO
          END DO
        END DO

        DO i=1,nx1-1
          DO j=1,ny1-1
            kk= MAX(1,MIN(lvlprof-1,                                    &
                INT((htrnx1y1(i,j)-zsnd(2))/snddelz)+2))
            IF(htrnx1y1(i,j) < zsnd(2)) k = 1
            tem = (htrnx1y1(i,j)-zsnd(kk))/(zsnd(kk+1)-zsnd(kk))
            ptbsfc(i,j)= ptsnd(kk)+(ptsnd(kk+1)-ptsnd(kk))*tem
            qvbsfc(i,j)= qvsnd(kk)+(qvsnd(kk+1)-qvsnd(kk))*tem
            pbsfc (i,j)= EXP(ALOG(psnd (kk))+                           &
                         (ALOG(psnd (kk+1))-ALOG(psnd (kk)))*tem)
          END DO
        END DO

      END IF

    END IF

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=ALOG(pbar(i,j,k)+pprt(i,j,k))
        END DO
      END DO
    END DO

    CALL intrpxyz3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                 &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    pprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,         &
                    temx1yz,temx1y1z)

    DO k=1,nz1-1
      DO j=1,ny1-1
        DO i=1,nx1-1
          pprt1(i,j,k)=EXP(pprt1(i,j,k))-pbar1(i,j,k)
        END DO
      END DO
    END DO
    DO j=1,ny1-1
      DO i=1,nx1-1
        ppsfc(i,j)=EXP(temx1y1z(i,j,2))-pbsfc(i,j)
      END DO
    END DO

    IF ( redo_base_state == 1 ) THEN

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem2(i,j,k)=ptprt(i,j,k)+ptbar(i,j,k)
          END DO
        END DO
      END DO

      CALL intrpxyz3d(tem2,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,               &
                      wgtsx,isx,wgtsy,jsy,wgtz,kz,                      &
                      intrphopt,intrpvopt,                              &
                      ptprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,      &
                      temx1yz,temx1y1z)

      DO j=1,ny1-1
        DO i=1,nx1-1
          ptpsfc(i,j)=temx1y1z(i,j,2)-ptbsfc(i,j)
        END DO
      END DO

      DO k=1,nz1-1
        DO j=1,ny1-1
          DO i=1,nx1-1
            ptprt1(i,j,k)=ptprt1(i,j,k)-ptbar1(i,j,k)
          END DO
        END DO
      END DO

    ELSE

      CALL intrpxyz3d(ptprt,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,              &
                     wgtsx,isx,wgtsy,jsy,wgtz,kz,                       &
                     intrphopt,intrpvopt,                               &
                     ptprt1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,       &
                     temx1yz,temx1y1z)

      DO j=1,ny1-1
        DO i=1,nx1-1
          ptpsfc(i,j)=temx1y1z(i,j,2)
        END DO
      END DO

    END IF

    CALL intrpxyz3d(qv,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qv1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qvsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(qc,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qc1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qcsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(qr,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qr1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qrsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(qi,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qi1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qisfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(qs,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qs1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qssfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(qh,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                   &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    qh1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,           &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        qhsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(tke,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                  &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    tke1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,          &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        tkesfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(kmh,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                  &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    kmh1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,          &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        kmhsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(kmv,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,                  &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    kmv1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,          &
                    temx1yz,temx1y1z)
    DO j=1,ny1-1
      DO i=1,nx1-1
        kmvsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO

    CALL intrpxyz3d(radfrc,nx,1,nx-1,ny,1,ny-1,nz,1,nz-1,               &
                    wgtsx,isx,wgtsy,jsy,wgtz,kz,                        &
                    intrphopt,intrpvopt,                                &
                    radfrc1,nx1,1,nx1-1, ny1,1,ny1-1,nz1,1,nz1-1,       &
                    temx1yz,temx1y1z)

    DO j=1,ny1-1
      DO i=1,nx1-1
        radsfc(i,j)=temx1y1z(i,j,2)
      END DO
    END DO
!
!  Reset extrapolated values below ground for certain options
!
    IF( npoint_below_ground /= 0 ) THEN ! Set values below input grid ground

      IF( bglopt == 2 ) THEN ! Constant extension below ground

        DO i=1,nx1-1
          DO j=1,ny1-1
            DO k=1,ktrnx1y1(i,j) ! only below ground
              ptprt1(i,j,k) = ptpsfc(i,j)
              pprt1 (i,j,k) = ppsfc (i,j)
              qv1   (i,j,k) = qvsfc (i,j)
              qc1   (i,j,k) = qcsfc (i,j)
              qr1   (i,j,k) = qrsfc (i,j)
              qi1   (i,j,k) = qisfc (i,j)
              qs1   (i,j,k) = qssfc (i,j)
              qh1   (i,j,k) = qhsfc (i,j)
              tke1  (i,j,k) = tkesfc(i,j)
              kmh1  (i,j,k) = kmhsfc(i,j)
              kmv1  (i,j,k) = kmhsfc(i,j)

              qvbar1(i,j,k) = qvbsfc(i,j)
              ptbar1(i,j,k) = ptbsfc(i,j)
              pbar1 (i,j,k) = pbsfc (i,j)
              rhobar1(i,j,k)= pbsfc(i,j)/                               &
                             (rd*ptbsfc(i,j)*(pbsfc(i,j)/p0)**rddcp)
              radfrc1(i,j,k)= radsfc(i,j)
            END DO
          END DO
        END DO

      ELSE IF( bglopt == 3 ) THEN ! Set to a missing value

        DO i=1,nx1-1
          DO j=1,ny1-1
            DO k=1,ktrnx1y1(i,j) ! only below ground
              ptprt1(i,j,k) = misvalue
              pprt1 (i,j,k) = misvalue
              qvbar1(i,j,k) = 0.0
              qv1   (i,j,k) = misvalue
              qc1   (i,j,k) = misvalue
              qr1   (i,j,k) = misvalue
              qi1   (i,j,k) = misvalue
              qs1   (i,j,k) = misvalue
              qh1   (i,j,k) = misvalue
              tke1  (i,j,k) = misvalue
              kmh1  (i,j,k) = misvalue
              kmv1  (i,j,k) = misvalue
              ptbar1(i,j,k) = 0.0
              pbar1 (i,j,k) = 0.0
              rhobar1(i,j,k)= misvalue
              radfrc1(i,j,k)= misvalue
            END DO
          END DO
        END DO
!
      ELSE IF( bglopt == 4  .OR. bglopt == 5) THEN 
         ! Using a constant lapse rate for temperature
         ! and then the hydrostatic balance for pressure

        DO i=1,nx1-1
          DO j=1,ny1-1
            DO k=1,ktrnx1y1(i,j) ! only below ground
              qv1   (i,j,k) = qvsfc (i,j)
              qc1   (i,j,k) = qcsfc (i,j)
              qr1   (i,j,k) = qrsfc (i,j)
              qi1   (i,j,k) = qisfc (i,j)
              qs1   (i,j,k) = qssfc (i,j)
              qh1   (i,j,k) = qhsfc (i,j)
              tke1  (i,j,k) = tkesfc(i,j)
              kmh1  (i,j,k) = kmhsfc(i,j)
              kmv1  (i,j,k) = kmhsfc(i,j)
              radfrc1(i,j,k)= radsfc(i,j)

!        ptprt1(i,j,k) = ptpsfc(i,j)
!        pprt1 (i,j,k) = ppsfc (i,j)

              qvbar1(i,j,k) = qvbsfc(i,j)

            END DO
          END DO
        END DO

        IF( redo_base_state == 1 ) THEN ! There are points below ground

          dtdz = -0.0068  ! Minus lapse rate

          DO i=1,nx1-1
            DO j=1,ny1-1

              k=ktrnx1y1(i,j)  ! First level below input grid ground
              IF( k >= 1) THEN
                delz = tem21(i,j,k)-htrnx1y1(i,j)
                ttotk1 = (ptbsfc(i,j)+ptpsfc(i,j)) *                    &
                        ((pbsfc (i,j)+ppsfc (i,j))/p0)**rddcp
                ttotk  = ttotk1+dtdz*delz
                lnptot = ALOG(pbsfc (i,j)+ppsfc (i,j))                  &
                         -g*2.0*delz/(rd*(ttotk1+ttotk))
                pprt1(i,j,k) = EXP(lnptot)-pbar1(i,j,k)
                ptprt1(i,j,k)= ttotk*(p0/(pbar1(i,j,k)+pprt1(i,j,k)))   &
                               **rddcp-ptbar1(i,j,k)
              END IF

              DO k=ktrnx1y1(i,j)-1,1,-1 ! other below-ground levels
                delz = tem21(i,j,k)-tem21(i,j,k+1)
                ttotk1 = (ptbar1(i,j,k+1)+ptprt1(i,j,k+1)) *            &
                        ((pbar1 (i,j,k+1)+pprt1 (i,j,k+1))/p0)**rddcp
                ttotk  = ttotk1+dtdz*delz
                lnptot = ALOG(pbar1(i,j,k+1)+pprt1(i,j,k+1))            &
                         -g*2.0*delz/(rd*(ttotk1+ttotk))
                pprt1(i,j,k) = EXP(lnptot)-pbar1(i,j,k)
                ptprt1(i,j,k)= ttotk*(p0/(pbar1(i,j,k)+pprt1(i,j,k)))   &
                               **rddcp-ptbar1(i,j,k)
              END DO
            END DO
          END DO

        END IF

      END IF

    END IF ! Check on npoint_below_ground

    800   CONTINUE

    xgrdorg = xorig1
    ygrdorg = yorig1

    CALL xytoll(1,1,xgrdorg,ygrdorg,origlat,origlon)
    CALL setorig( 2, origlat, origlon) ! set up the model origin on new grid

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
!    IF( nfile == 1) THEN
    IF(nfile == 1 .or. noutgrds > 0) THEN

      !wdt update
      DO i=1,nx1
        x1_out(i) = x1(i)-xgrdorg
      END DO

      !wdt update
      DO j=1,ny1
        y1_out(j) = y1(j)-ygrdorg
      END DO

      WRITE (*,*) "ctrlat & ctrlon:",ctrlat1,ctrlon1

      !wdt update: code block moved here
      !wdt begin block
      WRITE(6,'(/2x,a)')'For the output grid:'

      WRITE(6,'(2x,a,f15.7)')                                               &
          'The latitude  of the output grid center, ctrlat1=',ctrlat1
      WRITE(6,'(2x,a,f15.7/)')                                              &
          'The longitude of the output grid center, ctrlon1=',ctrlon1
      WRITE(6,'(2x,a,f15.7)')                                               &
          'The x coordinate of the output grid center, xctr1=',xctr1
      WRITE(6,'(2x,a,f15.7)')                                               &
          'The y coordinate of the output grid center, yctr1=',yctr1

      IF( ABS(xctr1-xctr1_old) > 0.0001*dx .OR.                             &
            ABS(yctr1-yctr1_old) > 0.0001*dy ) THEN
        WRITE(6,'(/1x,a/)')                                                 &
            '################################################################'
        WRITE(6,'(/2x,a/2x,a)')                                             &
            'Note that xctr1 and/or yctr1 had been changed by the program', &
            'due to snap_to_grid option. Check earlier messages.'
        WRITE(6,'(/1x,a/)')                                                 &
            '################################################################'
      END IF

      WRITE(6,'(/2x,a/2x,a,2f15.7,a)')                                      &
          'The SW corner (i,j)=(2,2) of the grid is located at ',           &
          '(',xgrdorg,ygrdorg,') of the input grid.'
      !wdt begin block

    END IF


!-----------------------------------------------------------------------
!
!  Copy grid information to the common block variables so that output
!  routines have the correct values for the new grid.
!
!-----------------------------------------------------------------------

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
    ! save global setting before setting them to new grid values
    ctrlat_sv = ctrlat
    ctrlon_sv = ctrlon
    latitud_sv = latitud
    dx_sv = dx
    dy_sv = dy
    dz_sv = dz
    dzmin_sv = dzmin
    strhopt_sv = strhopt
    zrefsfc_sv = zrefsfc
    dlayer1_sv = dlayer1
    dlayer2_sv = dlayer2
    zflat_sv = zflat
    strhtune_sv = strhtune
    nstyp_sv = nstyp

    ctrlat = ctrlat1
    ctrlon = ctrlon1
    latitud = ctrlat

    dx = dx1
    dy = dy1
    dz = dz1

    dzmin = dzmin1
    strhopt = strhopt1
    zrefsfc = zrefsfc1
    dlayer1 = dlayer11
    dlayer2 = dlayer21
    zflat = zflat1
    strhtune = strhtune1

    nstyp = nstyp1

    WRITE (cmnt(nocmnt),'(a,i4,a,i4,a,i4)')                             &
        ' nx =',nx1,', ny =',ny1,', nz =',nz1
!
!-----------------------------------------------------------------------
!
!  Print out the max/min of output varaibles.
!
!-----------------------------------------------------------------------
!
    WRITE(6,'(/1x,a/)')                                                 &
        'Min. and max. of input data interpolated to the new grid:'

    !wdt update
    CALL a3dmax0(x1_out,1,nx1,1,nx1,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(/1x,2(a,e13.6))') 'xmin    = ', amin,',  xmax    =',amax

    !wdt update
    CALL a3dmax0(y1_out,1,ny1,1,ny1,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'ymin    = ', amin,',  ymax    =',amax

    CALL a3dmax0(z1,1,nz1,1,nz1,1,1,1,1,1,1,1,1, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'zmin    = ', amin,',  zmax    =',amax

    CALL a3dmax0(zp1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1,           &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'zpmin   = ', amin,', zpmax    =',amax

    CALL a3dmax0(ubar1,1,nx1,1,nx1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,',  ubarmax =',amax

    CALL a3dmax0(vbar1,1,nx1,1,nx1-1,1,ny1,1,ny1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,',  vbarmax =',amax

    CALL a3dmax0(ptbar1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,      &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,',  ptbarmax=',amax

    CALL a3dmax0(pbar1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,       &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,',  pbarmax =',amax

    CALL a3dmax0(rhobar1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,     &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'rhobarmin=', amin,', rhobarmax=',amax

    CALL a3dmax0(qvbar1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,      &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,',  qvbarmax=',amax

    CALL a3dmax0(u1,1,nx1,1,nx1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,            &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'umin    = ', amin,',  umax    =',amax

    CALL a3dmax0(v1,1,nx1,1,nx1-1,1,ny1,1,ny1,1,nz1,1,nz1-1,            &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'vmin    = ', amin,',  vmax    =',amax

    CALL a3dmax0(w1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1,            &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'wmin    = ', amin,',  wmax    =',amax

    CALL a3dmax0(ptprt1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,      &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,',  ptprtmax=',amax

    CALL a3dmax0(pprt1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,       &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,',  pprtmax =',amax

    CALL a3dmax0(qv1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qvmin   = ', amin,',  qvmax   =',amax

    CALL a3dmax0(qc1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qcmin   = ', amin,',  qcmax   =',amax

    CALL a3dmax0(qr1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qrmin   = ', amin,',  qrmax   =',amax

    CALL a3dmax0(qi1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qimin   = ', amin,',  qimax   =',amax

    CALL a3dmax0(qs1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qsmin   = ', amin,',  qsmax   =',amax

    CALL a3dmax0(qh1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,         &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qhmin   = ', amin,',  qhmax   =',amax

    CALL a3dmax0(tke1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,        &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'tkemin  = ', amin,',  tkemax  =',amax

    CALL a3dmax0(kmh1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,        &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'kmhmin  = ', amin,',  kmhmax  =',amax

    CALL a3dmax0(kmv1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,        &
                amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'kmvmin  = ', amin,',  kmvmax  =',amax

    CALL a3dmax0(raing1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'raingmin= ', amin,',  raingmax=',amax

    CALL a3dmax0(rainc1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'raincmin= ', amin,',  raincmax=',amax

    CALL a3dmax0(prcrate1(1,1,1),1,nx1,1,nx1-1,1,ny1,1,ny1-1,           &
                 1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'prcr1min= ', amin,',  prcr1max=',amax

    CALL a3dmax0(prcrate1(1,1,2),1,nx1,1,nx1-1,1,ny1,1,ny1-1,           &
                 1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'prcr2min= ', amin,',  prcr2max=',amax

    CALL a3dmax0(prcrate1(1,1,3),1,nx1,1,nx1-1,1,ny1,1,ny1-1,           &
                 1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'prcr3min= ', amin,',  prcr3max=',amax

    CALL a3dmax0(prcrate1(1,1,4),1,nx1,1,nx1-1,1,ny1,1,ny1-1,           &
                 1,1,1,1,amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'prcr4min= ', amin,',  prcr4max=',amax

    DO iss = 0, nstyp1

      CALL a3dmax0(tsfc1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,  &
                   amax,amin)
      WRITE(6,'(1x,2(a,e13.6),a,i3)')                                   &
          'tsfcmin = ', amin,',  tsfcmax =',amax,' for soil type=',iss

      CALL a3dmax0(tsoil1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
                   amax,amin)
      WRITE(6,'(1x,2(a,e13.6),a,i3)')                                   &
          'tsoilmin= ', amin,',  tsoilmax=',amax,' for soil type=',iss

      CALL a3dmax0(wetsfc1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
                   amax,amin)
      WRITE(6,'(1x,2(a,e13.6),a,i3)')                                   &
          'wetsmin = ', amin,',  wetsmax =',amax,' for soil type=',iss

      CALL a3dmax0(wetdp1(1,1,iss),1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1, &
                   amax,amin)
      WRITE(6,'(1x,2(a,e13.6),a,i3)')                                   &
          'wetdmin = ', amin,',  wetdmax =',amax,' for soil type=',iss

      CALL a3dmax0(wetcanp1(1,1,iss),                                   &
                   1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6),a,i3)')                                   &
          'wetcmin = ', amin,',  wetcmax =',amax,' for soil type=',iss

    END DO

    CALL a3dmax0(roufns1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,           &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'roufnmin= ', amin,', roufnmax =',amax

    CALL a3dmax0(veg1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,              &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'vegmin  = ', amin,',   vegmax =',amax

    CALL a3dmax0(radfrc1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,nz1,1,nz1-1,     &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'radfnmin= ', amin,', radfnmax =',amax

    CALL a3dmax0(radsw1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,            &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'radswmin= ', amin,', radswmax =',amax

    CALL a3dmax0(rnflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,            &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'rnflxmin= ', amin,', rnflxmax =',amax

    CALL a3dmax0(usflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,            &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'usflxmin= ', amin,', usflxmax =',amax

    CALL a3dmax0(vsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,            &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'vsflxmin= ', amin,', vsflxmax =',amax

    CALL a3dmax0(ptsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,           &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'ptflxmin= ', amin,', ptflxmax =',amax

    CALL a3dmax0(qvsflx1,1,nx1,1,nx1-1,1,ny1,1,ny1-1,1,1,1,1,           &
                 amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qvflxmin= ', amin,', qvflxmax =',amax

    tem11= 0.0        ! To be put in place of wbar
!
!-----------------------------------------------------------------------
!
!  Data dump of the model grid and base state arrays:
!
!  First find a unique name basdmpfn(1:lbasdmpf) for the grid and
!  base state array dump file
!
!  If grid/base state data has been written out once, skip
!  the following writing block. Also no need to write out
!  separate data for Savi3D dump. The same for GrADS dump.
!
!-----------------------------------------------------------------------
!
    runname = new_runname
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
    IF (noutgrds > 0) THEN
      runname = name_grd(ng)
    ENDIF

    CALL gtlfnkey(runname, lfnkey)

    IF(houtfmt /= 9 ) THEN

      !wdt update
      IF( nfile > 1 ) GO TO 300 ! If done already, skip this part.

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

      PRINT*
      PRINT*,'Output grid/base state file is ', basdmpfn(1:lbasdmpf)

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

      !wdt update
      CALL setcornerll( nx1,ny1,x1_out,y1_out )
!wdt update
      CALL dtadump(nx1,ny1,nz1,nstyp1,hdmpfmt,nchout,                   &
                   basdmpfn(1:lbasdmpf),grdbas,filcmprs,                &
                   u1,v1,w1,ptprt1,pprt1,qv1,qc1,qr1,qi1,qs1,qh1,       &
                   tke1,kmh1,kmv1,                                      &
                   ubar1,vbar1,tem11,ptbar1,pbar1,rhobar1,qvbar1,       &
                   x1_out,y1_out,z1,zp1,hterain1,j11,j21,j31,  &
                   soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1,        &
                   tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1,      &
                   raing1,rainc1,prcrate1,                              &
                   radfrc1,radsw1,rnflx1,                               &
                   usflx1,vsflx1,ptsflx1,qvsflx1,                       &
                   tem21,wgtz,kz)

      !wdt update - removed gboutcnt

      300     CONTINUE

    END IF
!
!-----------------------------------------------------------------------
!
!  Then the time dependent fields:
!
!-----------------------------------------------------------------------
!
    IF( .NOT. (houtfmt == 9 .AND. vroutcnt == 1) ) THEN
!
!-----------------------------------------------------------------------
!
!  Reconstruct the file name using the specified directory name
!
!-----------------------------------------------------------------------
!
      CALL gtdmpfn(runname(1:lfnkey),dirname,                           &
                   ldirnam,curtim,hdmpfmt,1,0, hdmpfn, ldmpf)

    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
    IF ( exbcdmp == 4 ) rayklow = -1

    WRITE(6,'(a,a)') 'Writing t-dependent variable history dump ',      &
                      hdmpfn(1:ldmpf)
    grdbas = 0

!wdt update
    CALL dtadump(nx1,ny1,nz1,nstyp1,hdmpfmt,nchout,                     &
                 hdmpfn(1:ldmpf),grdbas,filcmprs,                       &
                 u1,v1,w1,ptprt1,pprt1,qv1,qc1,qr1,qi1,qs1,qh1,         &
                 tke1,kmh1,kmv1,                                        &
                 ubar1,vbar1,tem11,ptbar1,pbar1,rhobar1,qvbar1,         &
                 x1_out,y1_out,z1,zp1,hterain1,j11,j21,j31,  &
                 soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1,          &
                 tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1,        &
                 raing1,rainc1,prcrate1,                                &
                 radfrc1,radsw1,rnflx1,                                 &
                 usflx1,vsflx1,ptsflx1,qvsflx1,                         &
                 tem21,wgtz,kz)
!
!-----------------------------------------------------------------------
!
!  Write out soil model variable file
!
!-----------------------------------------------------------------------
!
    IF (sfcin == 1 .AND. soildmp > 0) THEN

      CALL cvttsnd( curtim, timsnd, tmstrln )

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

      IF( dirname /= ' ' ) THEN

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

      END IF

      CALL fnversn(soiloutfl, lfn)

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

      CALL wrtsoil(nx1,ny1,nstyp1, soiloutfl(1:lfn), dx1,dy1,           &
                   mapproj,trulat1,trulat2,trulon,sclfct,               &
                   ctrlat1,ctrlon1,                                     &
                   1,1,1,1,1,1,                                         &
                   tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,snowdpth1,soiltyp1)

      !wdt update
      IF (soildmp == 1) CALL soilcntl(nx1,ny1, soiloutfl(1:lfn),        &
                  1,1,1,1,1,1, x1_out,y1_out)

    END IF       ! sfcin.eq.1

    IF(nfile == 1) THEN  ! Need to do this only once
!
!-----------------------------------------------------------------------
!
!  Write out terrain data
!
!-----------------------------------------------------------------------
!
      IF (terndmp > 0) THEN
        CALL getunit( nunit )

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

        IF( dirname /= ' ' ) THEN

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

        END IF

        CALL fnversn(ternfn, lternfn )

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

        CALL writtrn(nx1,ny1,ternfn(1:lternfn), dx1,dy1,                &
                   mapproj,trulat1,trulat2,trulon,sclfct,               &
                   ctrlat1,ctrlon1,hterain1)

        !wdt update
        IF (terndmp == 1)                                               &
            CALL trncntl(nx1,ny1,ternfn(1:lternfn), x1_out,y1_out)

      END IF

!-----------------------------------------------------------------------
!
!  Write out surface property data file: sfcoutfl .
!
!-----------------------------------------------------------------------
!
      IF (landin == 1 .AND. sfcdmp > 0) THEN

        sfcoutfl = runname(1:lfnkey)//".sfcdata"
        lfn = lfnkey + 8

        IF( dirname /= ' ' ) THEN

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

        END IF

        CALL fnversn(sfcoutfl, lfn)

        PRINT *, 'Write surface property data in ',sfcoutfl(1:lfn)

        CALL wrtsfcdt(nx1,ny1,nstyp1,sfcoutfl(1:lfn), dx1,dy1,          &
                    mapproj,trulat1,trulat2,trulon,sclfct,              &
                    ctrlat1,ctrlon1,                                    &
                    1,1,1,1,1,0,                                        &
                    soiltyp1,stypfrct1,vegtyp1,lai1,roufns1,veg1,veg1)

        !wdt udpate
        IF (sfcdmp == 1) CALL sfccntl(nx1,ny1, sfcoutfl(1:lfn),         &
                     1,1,1,1,1,0, x1_out,y1_out)
      END IF       ! landin.eq.1

    END IF

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
    ctrlat = ctrlat_sv
    ctrlon = ctrlon_sv
    latitud = latitud_sv
    dx = dx_sv
    dy = dy_sv
    dz = dz_sv
    dzmin = dzmin_sv
    strhopt = strhopt_sv
    zrefsfc = zrefsfc_sv
    dlayer1 = dlayer1_sv
    dlayer2 = dlayer2_sv
    zflat = zflat_sv
    strhtune = strhtune_sv
    nstyp = nstyp_sv

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!wdt 2001-10-26 GMB multiple output grids
600 CONTINUE  !wdt
    IF (noutgrds > 0 .and. ng < noutgrds) GOTO 500

    vroutcnt = 1 ! Variables have been written out at least once

    IF (hinfmt == 9) GO TO 102

  END DO

  9001  CONTINUE

  !wdt update: moved code block to above

  WRITE (6,'(/a,g13.3,a)') "Current memory allocation=",                &
               current_memory_use,' (Mw).'
  WRITE (6,'(a,g13.3,a/)') "Maxumum memory allocation=",                &
               max_memory_use,' (Mw).'

  cpu_time = f_cputime() - cpu0
  PRINT*,'Total CPU time used =', cpu_time,' (s).'

  STOP

  9002  CONTINUE
  WRITE(6,'(1x,a,i2,/1x,a)')                                            &
      'Data read was unsuccessful. ireturn =', ireturn,                 &
      'Job stopped in ARPSINTRP.'

  STOP 9002
END PROGRAM arpsintrp


SUBROUTINE intrpx3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgtx,ix,          & 1
           intrphopt,aout,nx1,is1,ie1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Perform interpolation in the first dimension
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/1/1999.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
  INTEGER :: nx1,is1,ie1
  REAL :: ain (nx ,ny,nz)
  REAL :: aout(nx1,ny,nz)
  REAL :: wgtx(nx1,3)
  INTEGER :: ix(nx1)
  INTEGER :: intrphopt
  INTEGER :: i,j,k

  IF(intrphopt == 1) THEN
    DO k=ks ,ke
      DO j=js ,je
        DO i=is1,ie1
          aout(i,j,k)=      wgtx(i,1) *ain(ix(i)  ,j,k)                 &
                      +(1.0-wgtx(i,1))*ain(ix(i)+1,j,k)
        END DO
      END DO
    END DO
  ELSE
    DO k=ks ,ke
      DO j=js ,je
        DO i=is1,ie1
          aout(i,j,k)=wgtx(i,1)*ain(ix(i)-1,j,k)                        &
                     +wgtx(i,2)*ain(ix(i)  ,j,k)                        &
                     +wgtx(i,3)*ain(ix(i)+1,j,k)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE intrpx3d


SUBROUTINE intrpy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgty,jy,          & 1
           intrphopt,aout,ny1,js1,je1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Perform interpolation in the second dimension
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/1/1999.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
  INTEGER :: ny1,js1,je1
  REAL :: ain (nx,ny ,nz)
  REAL :: aout(nx,ny1,nz)
  REAL :: wgty(ny1,3)
  INTEGER :: jy(ny1)
  INTEGER :: intrphopt
  INTEGER :: i,j,k

  IF(intrphopt == 1) THEN
    DO k=ks ,ke
      DO j=js1,je1
        DO i=is ,ie
          aout(i,j,k)=      wgty(j,1) *ain(i,jy(j)  ,k)                 &
                      +(1.0-wgty(j,1))*ain(i,jy(j)+1,k)
        END DO
      END DO
    END DO
  ELSE
    DO k=ks ,ke
      DO j=js1,je1
        DO i=is ,ie
          aout(i,j,k)=wgty(j,1)*ain(i,jy(j)-1,k)                        &
                     +wgty(j,2)*ain(i,jy(j)  ,k)                        &
                     +wgty(j,3)*ain(i,jy(j)+1,k)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE intrpy3d


SUBROUTINE intrpz3d(ain,nx,is,ie, ny,js,je, nz,ks,ke, wgtz,kz,          & 1
           intrpvopt,aout,nz1,ks1,ke1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Perform interpolation in the third dimension
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/1/1999.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
  INTEGER :: nz1,ks1,ke1
  REAL :: ain (nx,ny,nz)
  REAL :: aout(nx,ny,nz1)
  REAL :: wgtz(nx,ny,nz1,3)
  INTEGER :: kz(nx,ny,nz1)
  INTEGER :: intrpvopt
  INTEGER :: i,j,k

  IF(intrpvopt == 1) THEN
    DO k=ks1,ke1
      DO i=is ,ie
        DO j=js ,je
          aout(i,j,k)=      wgtz(i,j,k,1) *ain(i,j,kz(i,j,k)  )         &
                      +(1.0-wgtz(i,j,k,1))*ain(i,j,kz(i,j,k)+1)
        END DO
      END DO
    END DO
  ELSE
    DO k=ks1,ke1
      DO i=is ,ie
        DO j=js ,je
          aout(i,j,k)=wgtz(i,j,k,1)*ain(i,j,kz(i,j,k)-1)                &
                     +wgtz(i,j,k,2)*ain(i,j,kz(i,j,k)  )                &
                     +wgtz(i,j,k,3)*ain(i,j,kz(i,j,k)+1)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE intrpz3d


SUBROUTINE intrpxy3d(ain,nx,is,ie, ny,js,je, nz,ks,ke,                  & 20,2
           wgtx,ix,wgty,jy,intrphopt,                                   &
           aout,nx1,is1,ie1, ny1,js1,je1,                               &
           temx1yz)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Perform interpolation in the first and second dimensions
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/1/1999.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
  INTEGER :: nx1,is1,ie1, ny1,js1,je1
  REAL :: ain (nx ,ny ,nz)
  REAL :: aout(nx1,ny1,nz)
  REAL :: wgtx(nx1,3),wgty(ny1,3)
  INTEGER :: ix(nx1),jy(ny1)
  INTEGER :: intrphopt

  REAL :: temx1yz(nx1,ny,nz)

  CALL intrpx3d(ain,nx,is,ie, ny,js,je, nz,ks,ke,                       &
                wgtx,ix,intrphopt, temx1yz,nx1,is1,ie1)

  CALL intrpy3d(temx1yz,nx1,is1,ie1, ny,js,je, nz,ks,ke,                &
                wgty,jy,intrphopt, aout,   ny1,js1,je1)

  RETURN
END SUBROUTINE intrpxy3d


SUBROUTINE intrpxyz3d(ain,nx,is,ie, ny,js,je, nz,ks,ke,                 & 22,2
           wgtx,ix,wgty,jy,wgtz,kz,                                     &
           intrphopt,intrpvopt,                                         &
           aout,nx1,is1,ie1, ny1,js1,je1, nz1,ks1,ke1,                  &
           temx1yz,temx1y1z)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Perform interpolation in all three dimensions
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/1/1999.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,is,ie, ny,js,je, nz,ks,ke
  INTEGER :: nx1,is1,ie1, ny1,js1,je1, nz1,ks1,ke1
  REAL :: ain (nx ,ny ,nz)
  REAL :: aout(nx1,ny1,nz1)
  REAL :: wgtx(nx1,3),wgty(ny1,3),wgtz(nx1,ny1,nz1,3)
  INTEGER :: ix(nx1),jy(ny1),kz(nx1,ny1,nz1)
  INTEGER :: intrphopt
  INTEGER :: intrpvopt

  REAL :: temx1yz (nx1,ny ,nz)
  REAL :: temx1y1z(nx1,ny1,nz)

  CALL intrpxy3d(ain, nx,is,ie, ny,js,je, nz,ks,ke,                     &
                 wgtx,ix,wgty,jy,intrphopt,                             &
                 temx1y1z, nx1,is1,ie1, ny1,js1,je1, temx1yz)

  CALL intrpz3d (temx1y1z, nx1,is1,ie1, ny1,js1,je1, nz,ks,ke,          &
                 wgtz,kz,intrpvopt,                                     &
                 aout, nz1,ks1,ke1)

  RETURN
END SUBROUTINE intrpxyz3d