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


subroutine heryin(pf      ,fint    ,fyb     ,fyt     ,y       , & 1,2
                  dy      ,ydp     ,jdp     ,fdp     ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! 
! Method: 
! For each departure point in the latitude slice to be forecast,
! interpolate (using unequally spaced Hermite cubic formulas) the
! x interpolants to the y value of the departure point.
! 
! Author:
! Original version:  J. Olson
! Standardized:      J. Rosinski, June 1992
! Reviewed:          D. Williamson, P. Rasch, August 1992
! Reviewed:          D. Williamson, P. Rasch, March 1996
!
!-----------------------------------------------------------------------
!
! $Id: heryin.F90,v 1.1.44.1 2002/06/15 13:47:45 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
!-----------------------------------------------------------------------
   implicit none
!------------------------------Parameters-------------------------------
#include <parslt.h>
!------------------------------Arguments--------------------------------
!
! Input arguments
!
   integer, intent(in) :: pf                     ! dimension (number of fields)
!
   real(r8), intent(in) :: fint(plon,plev,ppdy,pf)   ! x-interpolants
   real(r8), intent(in) :: fyb (plon,plev,pf)        ! y-derivatives at bottom of interval
   real(r8), intent(in) :: fyt (plon,plev,pf)        ! y-derivatives at top    of interval
   real(r8), intent(in) :: y   (platd)               ! latitude grid coordinates
   real(r8), intent(in) :: dy  (platd)               ! intervals between latitude grid pts.
   real(r8), intent(in) :: ydp (plon,plev)           ! lat. coord of departure point.
!
   integer, intent(in) :: jdp (plon,plev)        ! lat. index of departure point.
   integer, intent(in) :: nlon
!
! Output arguments
!
   real(r8), intent(out) :: fdp (plon,plev,pf)        ! y-interpolants

!
!-----------------------------------------------------------------------
!
!  pf      Number of fields being interpolated.
!  fint    (fint(i,k,j,m),j=ppdy/2,ppdy/2 + 1) contains the x
!          interpolants at the endpoints of the y-interval that
!          contains the departure point for grid point (i,k).  The last
!          index of fint allows for interpolation of multiple fields.
!          fint is generated by a call to herxin.
!  fyb     fyb(i,k,.) is the derivative at the "bottom" of the
!          y-interval that contains the departure point of grid
!          point (i,k).  fyb is generated by a call to cubydr.
!  fyt     fyt(i,k,.) is the derivative at the "top" of the y-interval
!          that contains the departure point of grid point (i,k).
!          fyt is generated by a call to cubydr.
!  y       y-coordinate (latitude) values in the extended array.
!  dy      Increment in the y-coordinate value for each interval in the
!          extended array.
!  ydp     ydp(i,k) is the y-coordinate of the departure point that
!          corresponds to global grid point (i,k) in the latitude slice
!          being forecasted.
!  jdp     jdp(i,k) is the index of the y-interval that contains the
!          departure point corresponding to global grid point (i,k) in
!          the latitude slice being forecasted.
!          Note that
!                y(jdp(i,k)) .le. ydp(i,k) .lt. y(jdp(i,k)+1) .
!  fdp     Horizontally interpolated field values at the departure point
!          for the latitude slice being forecasted.
!
!---------------------------Local variables-----------------------------
!
   integer i,k               ! index
   integer jb                ! index corresponding to bot of interval
   integer jt                ! index corresponding to top of interval
   integer m                 ! index
!
   real(r8) dyj(plon,plev)   ! latitude interval containing dep. pt.
   real(r8) yb (plon,plev)   ! |
   real(r8) yt (plon,plev)   ! |
   real(r8) hb (plon,plev)   ! | -- interpolation coefficients
   real(r8) ht (plon,plev)   ! |
   real(r8) dhb(plon,plev)   ! |
   real(r8) dht(plon,plev)   ! |
!
!-----------------------------------------------------------------------
!
   jb = ppdy/2
   jt = jb + 1
!
   do k=1,plev
      do i = 1,nlon
         dyj(i,k) = dy(jdp(i,k))
         yb (i,k) = ( y(jdp(i,k)+1) - ydp(i,k) )/dyj(i,k)
         yt (i,k) = 1. - yb(i,k)
         hb (i,k) = ( 3.0 - 2.0*yb(i,k) )*yb(i,k)**2
         ht (i,k) = ( 3.0 - 2.0*yt(i,k) )*yt(i,k)**2
         dhb(i,k) = -dyj(i,k)*( yb(i,k) - 1. )*yb(i,k)**2
         dht(i,k) =  dyj(i,k)*( yt(i,k) - 1. )*yt(i,k)**2
      end do
   end do
!
! Loop over fields.
!
   do m = 1,pf
      do k=1,plev
         do i = 1,nlon
            fdp(i,k,m) = fint(i,k,jb,m)*hb(i,k) + fyb(i,k,m)*dhb(i,k) + &
               fint(i,k,jt,m)*ht(i,k) + fyt(i,k,m)*dht(i,k)
         end do
      end do
   end do
!
   return
end subroutine heryin