#include <misc.h>


subroutine lininterp (arrin, yin, nlev, nlatin, arrout, & 6,5
                      yout, nlatout)
!----------------------------------------------------------------------- 
! 
! Purpose: Do a linear interpolation from input mesh defined by yin to output
!          mesh defined by yout.  Where extrapolation is necessary, values will
!          be copied from the extreme edge of the input grid.  Vectorization is over
!          the vertical (nlev) dimension.
! 
! Method: Check validity of input, then determine weights, then do the N-S interpolation.
! 
! Author: Jim Rosinski
!-----------------------------------------------------------------------
   use shr_kind_mod, only: r8 => shr_kind_r8
   use abortutils, only: endrun
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
!
! Arguments
!
   integer, intent(in) :: nlev                   ! number of vertical levels
   integer, intent(in) :: nlatin                 ! number of input latitudes
   integer, intent(in) :: nlatout                ! number of output latitudes

   real(r8), intent(in) :: arrin(nlev,nlatin)    ! input array of values to interpolate
   real(r8), intent(in) :: yin(nlatin)           ! input mesh
   real(r8), intent(in) :: yout(nlatout)         ! output mesh

   real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array
!
! Local workspace
!
   integer j, jj              ! latitude indices
   integer js, jn, jjprev     ! latitude indices
   integer k                  ! level index
   integer icount             ! number of values

   real(r8) extrap            ! percent grid non-overlap
!
! Dynamic
!
   integer :: jjm(nlatout)
   integer :: jjp(nlatout)

   real(r8) :: wgts(nlatout)
   real(r8) :: wgtn(nlatout)
!
! Check validity of input coordinate arrays: must be monotonically increasing,
! and have a total of at least 2 elements
!
   if (nlatin.lt.2) then
      call endrun ('LININTERP: Must have at least 2 input points for interpolation')
   end if

   icount = 0
   do j=1,nlatin-1
      if (yin(j).gt.yin(j+1)) icount = icount + 1
   end do

   do j=1,nlatout-1
      if (yout(j).gt.yout(j+1)) icount = icount + 1
   end do

   if (icount.gt.0) then
      call endrun ('LININTERP: Non-monotonic coordinate array(s) found')
   end if
!
! Initialize index arrays for later checking
!
   do j=1,nlatout
      jjm(j) = 0
      jjp(j) = 0
   end do
!
! For values which extend beyond N and S boundaries, set weights
! such that values will just be copied.
!
   do js=1,nlatout
      if (yout(js).gt.yin(1)) goto 10
      jjm(js) = 1
      jjp(js) = 1
      wgts(js) = 1.
      wgtn(js) = 0.
   end do
10 do jn=nlatout,1,-1
      if (yout(jn).le.yin(nlatin)) goto 20
      jjm(jn) = nlatin
      jjp(jn) = nlatin
      wgts(jn) = 1.
      wgtn(jn) = 0.
   end do
!
! Loop though output indices finding input indices and weights
!
20 jjprev = 1
   do j=js,jn
      do jj=jjprev,nlatin-1
         if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
            jjm(j) = jj
            jjp(j) = jj + 1
            wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
            wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
            goto 30
         end if
      end do
      write(6,*)'LININTERP: Failed to find interp values'
30    jjprev = jj
   end do
!
! Check grid overlap
!
   extrap = 100.*((js - 1.) + float(nlatout - jn))/nlatout
   if (extrap.gt.30.) then
      write(6,*)'LININTERP WARNING:',extrap,' % of output grid will have to be extrapolated'
   end if
!
! Check that interp/extrap points have been found for all outputs
!
   icount = 0
   do j=1,nlatout
      if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1
   end do
   if (icount.gt.0) then
      call endrun ('LININTERP: Point found without interp indices')
   end if
!
! Do the interpolation
!
   do j=1,nlatout
      do k=1,nlev
         arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j)
      end do
   end do

   return
end subroutine lininterp