subroutine limdy(pf      ,fint    ,dy      ,jdp     ,fyb     ,&
                 fyt     ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Limit the y-derivative estimates so they satisy the SCM0 for the
! x-interpolated data corresponding to the departure points of a single
! latitude slice in the global grid, that is, they are monotonic, but
! spline has only C0 continuity
! 
! 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$
! $Author$
!
!-----------------------------------------------------------------------

   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid,       only: plon, plev
   use scanslt,      only: platd
!-----------------------------------------------------------------------
   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) :: dy(platd)               ! interval lengths in lat grid
!
   integer, intent(in) :: jdp(plon,plev)       ! j-index of coord. of dep. pt.
   integer, intent(in) :: nlon
!
! Input/output arguments
!
   real(r8), intent(inout) :: fyb(plon,plev,pf)       ! y-derivatives at bot of interval
   real(r8), intent(inout) :: fyt(plon,plev,pf)       ! y-derivatives at top 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.
!  dy      Increment in the y-coordinate value for each interval in the
!          extended array.
!  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 limited derivative at the bot of the y
!          interval that contains the departure point of global grid
!          point (i,k).
!  fyt     fyt(i,k,.) is the limited derivative at the top of the y
!          interval that contains the departure point of global grid
!          point (i,k).
!
!---------------------------Local variables-----------------------------
!
   integer i,k,m             ! indices
   integer jb                ! index for bottom of interval
   integer jt                ! index for top    of interval
!
   real(r8) rdy (plon,plev)      ! 1./dy
   real(r8) deli(plon)           ! simple linear derivative

!GRCJR
   real(r8) fac,tmp1,tmp2
   fac = 3._r8*(1._r8 - 10._r8*epsilon(fac))
!
!-----------------------------------------------------------------------
!
   jb = ppdy/2
   jt = jb + 1
!
!$OMP PARALLEL DO PRIVATE (K, I)
   do k = 1,plev
      do i = 1,nlon
         rdy(i,k) = 1._r8/dy(jdp(i,k))
      end do
   end do
!
! Loop over fields.
!
   do m = 1,pf
!$OMP PARALLEL DO PRIVATE (K, I, DELI, TMP1, TMP2)
      do k = 1,plev
         do i = 1,nlon
            deli(i) = ( fint(i,k,jt,m) - fint(i,k,jb,m) )*rdy(i,k)
!         end do
!
! Limiter
!
!GRCJR         call scm0(nlon,deli,fyb(1,k,m),fyt(1,k,m))
!         do i = 1,nlon
             tmp1 = fac*deli(i)
             tmp2 = abs( tmp1 )
             if( deli(i)*fyb(i,k,m)   <= 0.0_r8  ) fyb(i,k,m) = 0._r8
             if( deli(i)*fyt(i,k,m)   <= 0.0_r8  ) fyt(i,k,m) = 0._r8
             if( abs(    fyb(i,k,m) ) >  tmp2 ) fyb(i,k,m) = tmp1
             if( abs(    fyt(i,k,m) ) >  tmp2 ) fyt(i,k,m) = tmp1
         end do
      end do
   end do
!
   return
end subroutine limdy