#include <misc.h>
#include <params.h>


subroutine vertinterp(ncol, ncold, nlev, pmid, pout, arrin, arrout) 18,3

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Vertically interpolate input array to output pressure level
! Copy values at boundaries.
! 
! Method: 
! 
! Author: 
! 
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use abortutils,   only: endrun

  implicit none

!------------------------------Arguments--------------------------------
  integer , intent(in)  :: ncol              ! column dimension
  integer , intent(in)  :: ncold             ! declared column dimension
  integer , intent(in)  :: nlev              ! vertical dimension
  real(r8), intent(in)  :: pmid(ncold,nlev)  ! input level pressure levels 
  real(r8), intent(in)  :: pout              ! output pressure level 
  real(r8), intent(in)  :: arrin(ncold,nlev) ! input  array
  real(r8), intent(out) :: arrout(ncold)     ! output array (interpolated)
!--------------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer i,k               ! indices
  integer kupper(ncold)     ! Level indices for interpolation
  real(r8) dpu              ! upper level pressure difference
  real(r8) dpl              ! lower level pressure difference
  logical found(ncold)      ! true if input levels found
  logical error             ! error flag 
!-----------------------------------------------------------------
!
! Initialize index array and logical flags
!
  do i=1,ncol
     found(i)  = .false.
     kupper(i) = 1
  end do
  error = .false.
!     
! Store level indices for interpolation. 
! If all indices for this level have been found, 
! do the interpolation 
!     
  do k=1,nlev-1
     do i=1,ncol
        if ((.not. found(i)) .and. pmid(i,k)<pout .and. pout<=pmid(i,k+1)) then
           found(i) = .true.
           kupper(i) = k
        end if
     end do
  end do
!
! If we've fallen through the k=1,nlev-1 loop, we cannot interpolate and
! must extrapolate from the bottom or top data level for at least some
! of the longitude points.
!
  do i=1,ncol
     if (pout <= pmid(i,1)) then
        arrout(i) = arrin(i,1)
     else if (pout >= pmid(i,nlev)) then
        arrout(i) = arrin(i,nlev)
     else if (found(i)) then
        dpu = pout - pmid(i,kupper(i))
        dpl = pmid(i,kupper(i)+1) - pout
        arrout(i) = (arrin(i,kupper(i)  )*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu)
     else
        error = .true.
     end if
  end do
!     
! Error check
!
  if (error) then
     call endrun ('VERTINTERP: ERROR FLAG')
  end if

  return
end subroutine vertinterp