#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