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


subroutine lagyin(pf      ,fint    ,wdy     ,ydp     ,jdp     ,  & 1,7
                  jcen    ,fdp     ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! For each departure point in the latitude slice to be forecast,
! interpolate (using unequally spaced Lagrange cubic formulas) the
! x interpolants to the y value of the departure point.
! 
! Method: 
! 
! 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: lagyin.F90,v 1.1.44.3 2004/04/28 23:43:55 eaton Exp $
! $Author: eaton $
!
!-----------------------------------------------------------------------

   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use abortutils, only: endrun
#if (!defined UNICOSMP)
   use srchutil, only: whenieq
#endif
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#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) :: wdy(4,2,platd)           ! y-interpolation weights
   real(r8), intent(in) :: ydp(plon,plev)           ! y-coordinates of departure pts.
!
   integer, intent(in) :: jdp(plon,plev)            ! j-index of departure point coord.
   integer, intent(in) :: jcen                      ! current latitude
   integer, intent(in) :: nlon
!
! Output arguments
!
   real(r8), intent(out) :: fdp(plon,plev,pf)         ! interpolants at the horiz. depart. pt.
!
!-----------------------------------------------------------------------
!
!  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.
!  wdy     Grid values and weights for Lagrange cubic interpolation on
!          the unequally spaced y-grid.
!  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,m               ! indices
!
   real(r8) ymy1                 ! |
   real(r8) ymy2                 ! |
   real(r8) ymy3                 ! |
   real(r8) ymy4                 ! |
   real(r8) coef12               ! |
   real(r8) coef34               ! | -- interpolation weights/coeffs.
   real(r8) term1(plon,plev)          ! |
   real(r8) term2(plon,plev)          ! |
   real(r8) term3(plon,plev)          ! |
   real(r8) term4(plon,plev)          ! |
!
   integer jdpval,icount,ii,indx(plon),nval
   integer k
!
!-----------------------------------------------------------------------
!
   if( ppdy .ne. 4) then
      call endrun ('LAGYIN:Error:  ppdy .ne. 4')
   end if
   icount = 0
   do jdpval=jcen-2,jcen+1
      if (icount.lt.nlon*plev) then
         do k=1,plev
            call whenieq(nlon,jdp(1,k),1,jdpval,indx,nval)
            icount = icount + nval
!
            do ii = 1,nval
               i=indx(ii)
               ymy3     = ydp(i,k) - wdy(3,1,jdpval)
               ymy4     = ydp(i,k) - wdy(4,1,jdpval)
               coef12   = ymy3*ymy4
               ymy2     = ydp(i,k) - wdy(2,1,jdpval)
               term1(i,k) = coef12*ymy2*wdy(1,2,jdpval)
               ymy1     = ydp(i,k) - wdy(1,1,jdpval)
               term2(i,k) = coef12*ymy1*wdy(2,2,jdpval)
               coef34   = ymy1*ymy2   
               term3(i,k) = coef34*ymy4*wdy(3,2,jdpval)
               term4(i,k) = coef34*ymy3*wdy(4,2,jdpval)
            end do
         end do
      end if
   end do
   if (icount.ne.nlon*plev) then
      write(6,*)'LAGYIN: Departure pt out of bounds: jcen,icount,nlon*plev=',jcen,icount,nlon*plev
      call endrun
   end if
!
! Loop over fields.
!
   do m = 1,pf
      do k=1,plev
         do i = 1,nlon
            fdp(i,k,m) = fint(i,k,1,m)*term1(i,k) +  &
               fint(i,k,2,m)*term2(i,k) +  &
               fint(i,k,3,m)*term3(i,k) +  &
               fint(i,k,4,m)*term4(i,k)
         end do
      end do
   end do
!
   return
end subroutine lagyin