#include "dims.h"
      SUBROUTINE DLDP(NX,XL,XP,LEN)
      use cons_module,only: dlamda_1div12,dlamda_2div3,
     |  dphi_2div3,dphi_1div12
      implicit none
#include "params.h"
#include "fcom.h"
#include "buff.h"
#include "phys.h"
!
! Args:
!     real,intent(out) :: XL(1),XP(1)
      real,intent(out) :: XL(zimxp,zkmxp),XP(zimxp,zkmxp)
      integer,intent(in) :: nx,len
!
! Local:
      integer :: lenp2,nxk,i,nxjp2k,nxjp1k,nxjm1k,nxjm2k
      integer :: k
!     real,dimension(zimxp,zkmxp) :: fj,
!    |  fip1,fim1,fip2,fim2, fjp1,fjm1,fjp2,fjm2
!
!     call addfsech('DLDP_ZIN',' ',' ',f(1,nj+nx),zimxp,zkmxp,zkmx,j)

      LENP2=LEN+2 ! equiv to len3-4+2 = len3-2, i.e., k=1,zkmxp, i=3,74
C     ****                   COMPUTE LONGITUDINAL DERIVATIVE
      NXK   = NJ+NX

!     DO 10 I=3,LENP2
!     XL(I)=dlamda_2div3*(F(I+1,NXK)-F(I-1,NXK))-dlamda_1div12*
!    |  (F(I+2,NXK)-F(I-2,NXK))
!  10 CONTINUE

      do k=1,zkmxp
        do i=3,zimxp-2
          XL(I,k)=dlamda_2div3 *(F(I+1,NXK)-F(I-1,NXK))-
     |            dlamda_1div12*(F(I+2,NXK)-F(I-2,NXK))

!         fip1(i,k) = f(i+1,nxk)
!         fim1(i,k) = f(i-1,nxk)
!         fip2(i,k) = f(i+2,nxk)
!         fim2(i,k) = f(i-2,nxk)

        enddo ! i=3,zimxp-2
        nxk = nxk+1
      enddo ! k=1,zkmxp

!     do i=3,zimxp-2
!       write(6,"('dldp zl: lat=',i2,' i=',i2)") j,i
!       write(6,"('f(i+1,:)=',/,(6e12.4))") fip1(i,:)
!       write(6,"('f(i-1,:)=',/,(6e12.4))") fim1(i,:)
!       write(6,"('f(i+2,:)=',/,(6e12.4))") fip2(i,:)
!       write(6,"('f(i-2,:)=',/,(6e12.4))") fim2(i,:)
!       write(6,"('zl(i,:)=',/,(6e12.4))") xl(i,:)
!     enddo


C     ****                   COMPUTE LATITUDINAL DERIVATIVE
      nxk = nj+nx
      NXJP2K= NJP2+NX
      NXJP1K= NJP1+NX
      NXJM1K= NJM1+NX
      NXJM2K= NJM2+NX

!     DO 20 I=3,LENP2
!     XP(I)=dphi_2div3*(F(I,NXJP1K)-F(I,NXJM1K))-
!    1      dphi_1div12*(F(I,NXJP2K)-F(I,NXJM2K))
!  20 CONTINUE

      do k=1,zkmxp
!       do i=1,zimxp
!         fj(i,k)   = f(i,nxk)
!         fjp1(i,k) = f(i,nxjp1k)
!         fjm1(i,k) = f(i,nxjm1k)
!         fjp2(i,k) = f(i,nxjp2k)
!         fjm2(i,k) = f(i,nxjm2k)
!       enddo
        nxk = nxk+1
        NXJP2K= nxjp2k+1
        NXJP1K= nxjp1k+1
        NXJM1K= nxjm1k+1
        NXJM2K= nxjm2k+1
      enddo

      nxk = nj+nx
      NXJP2K= NJP2+NX
      NXJP1K= NJP1+NX
      NXJM1K= NJM1+NX
      NXJM2K= NJM2+NX
      do k=1,zkmxp
        do i=3,zimxp-2

          XP(I,k)=dphi_2div3 *(F(I,NXJP1K)-F(I,NXJM1K))-
     1            dphi_1div12*(F(I,NXJP2K)-F(I,NXJM2K))
        enddo ! i=3,zimxp-2

        nxk = nxk+1
        NXJP2K= nxjp2k+1
        NXJP1K= nxjp1k+1
        NXJM1K= nxjm1k+1
        NXJM2K= nxjm2k+1
      enddo ! k=1,zkmxp

!     do i=3,zimxp-2
!       write(6,"('dldp zp: lat=',i2,' i=',i2,' zin=',/,(6e12.4))")
!    |    j,i,fj(i,:)
!       write(6,"('dldp zp: lat=',i2,' i=',i2)") j,i
!       write(6,"('zin(:,i,lat+1)=',/,(6e12.4))") fjp1(i,:)
!       write(6,"('zin(:,i,lat-1)=',/,(6e12.4))") fjm1(i,:)
!       write(6,"('zin(:,i,lat+2)=',/,(6e12.4))") fjp2(i,:)
!       write(6,"('zin(:,i,lat-2)=',/,(6e12.4))") fjm2(i,:)
!       write(6,"('zp(:,i,lat)=',/,(6e12.4))") xp(i,:)
!     enddo


!     call addfsech('Z_JM2'  ,' ',' ',fjm2,zimxp,zkmxp,zkmx,j)
!     call addfsech('Z_JM1'  ,' ',' ',fjm1,zimxp,zkmxp,zkmx,j)
!     call addfsech('Z_JP1'  ,' ',' ',fjp1,zimxp,zkmxp,zkmx,j)
!     call addfsech('Z_JP2'  ,' ',' ',fjp2,zimxp,zkmxp,zkmx,j)

!     call addfsech('DLDP_ZL',' ',' ',xl,zimxp,zkmxp,zkmx,j)
!     call addfsech('DLDP_ZP',' ',' ',xp,zimxp,zkmxp,zkmx,j)

      RETURN
      END
C