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

subroutine cubydr(pf      ,fint    ,wdy     ,jdp     ,jcen    , & 1,6
                  fyb     ,fyt     ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute Lagrangian cubic derivative estimates at both ends of the
! intervals in the y coordinate (unequally spaced) containing the
! departure points for the latitude slice being forecasted.
!
! 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: cubydr.F90,v 1.1.44.3 2004/04/28 23:43:53 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                     ! number of constituent fields
!
  real(r8), intent(in) :: fint(plon,plev,ppdy,pf)   ! constituent x- interpolants
  real(r8), intent(in) :: wdy(4,2,platd)            ! latitude interpolation weights
!
  integer, intent(in) :: jdp(plon,plev)         ! indices of latitude intervals
  integer, intent(in) :: jcen                   ! current latitude index
  integer, intent(in) :: nlon
!
! Output arguments
!
  real(r8), intent(out) :: fyb(plon,plev,pf)         ! Derivative at south end of interval
  real(r8), intent(out) :: fyt(plon,plev,pf)         ! Derivative at north end of interval
!-----------------------------------------------------------------------
!
!  pf      Number of fields being interpolated.
!  fint    (fint(i,k,j,m),j=1,ppdy) contains the x interpolants at each
!          latitude needed for the y derivative estimates at the
!          endpoints of the 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     Weights for Lagrange cubic derivative estimates on the
!          unequally spaced latitude grid.  If grid interval j (in
!          extended array) is surrounded by a 4 point stencil, then
!          the derivative at the "bottom" of the interval uses the
!          weights wdy(1,1,j),wdy(2,1,j), wdy(3,1,j), and wdy(4,1,j).
!          The derivative at the "top" of the interval uses wdy(1,2,j),
!          wdy(2,2,j), wdy(3,2,j), and wdy(4,2,j).
!  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.
!          Suppose yb contains the y-coordinates of the extended array
!          and ydp(i,k) is the y-coordinate of the departure point
!          corresponding to grid point (i,k).  Then,
!                yb(jdp(i,k)) .le. ydp(i,k) .lt. yb(jdp(i,k)+1) .
!  fyb     fyb(i,k,.) is the derivative at the bottom of the y interval
!          that contains the departure point of global grid point (i,k).
!  fyt     fyt(i,k,.) is the derivative at the top of the y interval
!          that contains the departure point of global grid point (i,k).
!
!---------------------------Local variables-----------------------------
!
  integer i,k               ! index
  integer m                 ! index
  integer jdpval            ! index
  integer icount            ! counter
  integer ii                ! index
  integer indx(plon)        ! set of indices for indirect addressing
  integer nval              ! number of indices for given "jdpval"
!
!-----------------------------------------------------------------------
!
  icount = 0
  do jdpval=jcen-2,jcen+1
     do k=1,plev
        call whenieq(nlon,jdp(1,k),1,jdpval,indx,nval)
        icount = icount + nval
!DIR$ PREFERSTREAM
        do m=1,pf
!DIR$ CONCURRENT
!DIR$ PREFERVECTOR
           do ii=1,nval
              i=indx(ii)
              fyb(i,k,m) = wdy(1,1,jdpval)*fint(i,k,1,m) +  &
                   wdy(2,1,jdpval)*fint(i,k,2,m) +  &
                   wdy(3,1,jdpval)*fint(i,k,3,m) +  &
                   wdy(4,1,jdpval)*fint(i,k,4,m)
!
              fyt(i,k,m) = wdy(1,2,jdpval)*fint(i,k,1,m) +  &
                   wdy(2,2,jdpval)*fint(i,k,2,m) +  &
                   wdy(3,2,jdpval)*fint(i,k,3,m) +  &
                   wdy(4,2,jdpval)*fint(i,k,4,m)
           end do
        end do
     end do
     if (icount.eq.nlon*plev) return
  end do
  if (icount.ne.nlon*plev) then
     write(6,*)'CUBYDR: Departure point out of bounds: jcen,icount,nlon*plev=',jcen,icount,nlon*plev
     call endrun ()
  end if
!
  return
end subroutine cubydr