!> @author Ignis
    !> @brief 특정한 order of accuracy를 가지는 compact finite difference scheme (generalized Padé scheme)을 구현한 모듈입니다.
    !!
    !! 이 모듈은 high-order compact 차분 스킴을 정의하고 수치 미분을 수행합니다.
    !! Compact 차분 스킴은 implicit 스킴이므로, 모듈 내부에서 implicit 연산을 빠르게 해결하기 위해 
    !! tridiagonal matrix solver(stdlu, ptdlu 등)가 함께 구현되어 동작합니다.
    MODULE Compact
     use, intrinsic :: iso_fortran_env, only: real64
     IMPLICIT NONE

     REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & !< x-방향 LU Decomposition 밴드 계수
                                            lyf,lys,wyf,wys, & !< y-방향 LU Decomposition 밴드 계수
                                            lzf,lzs,wzf,wzs    !< z-방향 LU Decomposition 밴드 계수

     INTEGER :: nxc,nyc,nzc !< 각 방향의 실제 격자 사이즈 수치 (x, y, z)

     REAL(KIND=8), PARAMETER :: ezero = 1.0e-14
      
    CONTAINS

     !> Entry point for LU decomposition calculations.
     !!
     !! Prepares the workspace allocations and runs the decomposition calculation for all three directions.
     !!
     !! @param nx Grid size in x-direction.
     !! @param ny Grid size in y-direction.
     !! @param nz Grid size in z-direction.
     !! @param xp Periodic flag for x-direction (0 = periodic, other = non-periodic).
     !! @param yp Periodic flag for y-direction (0 = periodic, other = non-periodic).
     !! @param zp Periodic flag for z-direction (0 = periodic, other = non-periodic).
     SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp)
      INTEGER, INTENT(IN) :: nx,ny,nz
      INTEGER, INTENT(IN) :: xp,yp,zp
      INTEGER :: ierr

       nxc=nx
       nyc=ny
       nzc=nz

       CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp)

       CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp)

     END SUBROUTINE ludcmp

     !> LU 분해를 위한 포트란 workspace 배열 메모리를 동적 할당하는 서브루틴입니다.
     !!
     !! 경계 조건(xp, yp, zp = 0 주기적 경계 조건, 1 비주기적 경계 조건)에 따라 배열 크기와 
     !! 할당 여부를 결정하며, 할당에 실패하면 에러를 출력하고 즉시 프로그램을 안전하게 종료(STOP)시킵니다.
     !!
     !! @param nx Grid size in x-direction.
     !! @param ny Grid size in y-direction.
     !! @param nz Grid size in z-direction.
     !! @param xp Periodic flag for x-direction (0 = periodic, other = non-periodic).
     !! @param yp Periodic flag for y-direction (0 = periodic, other = non-periodic).
     !! @param zp Periodic flag for z-direction (0 = periodic, other = non-periodic).
     SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp)
      INTEGER, INTENT(IN) :: nx,ny,nz
      INTEGER, INTENT(IN) :: xp,yp,zp
      INTEGER :: ierr

       nxc=nx
       nyc=ny
       nzc=nz

!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'

!      xp, yp, zp = 0 : periodic
       ALLOCATE(lxf(nxc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       ALLOCATE(lxs(nxc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       IF(xp.eq.0) THEN
         ALLOCATE(wxf(nxc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
         ALLOCATE(wxs(nxc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
       ENDIF

       ALLOCATE(lyf(nyc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       ALLOCATE(lys(nyc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       IF(yp.eq.0) THEN
         ALLOCATE(wyf(nyc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
         ALLOCATE(wys(nyc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
       ENDIF

       ALLOCATE(lzf(nzc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       ALLOCATE(lzs(nzc),STAT=ierr)
       IF(ierr /= 0) THEN
         PRINT*, 'work array for lud allocation failed'
         STOP 1
       ENDIF
       IF(zp.eq.0) THEN
         ALLOCATE(wzf(nzc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
         ALLOCATE(wzs(nzc),STAT=ierr)
         IF(ierr /= 0) THEN
           PRINT*, 'work array for lud allocation failed'
           STOP 1
         ENDIF
       ENDIF

     END SUBROUTINE ludcmp_allocate

     SUBROUTINE ludcmp_deallocate(xp,yp,zp)
      INTEGER, INTENT(IN) :: xp,yp,zp

!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'

!      xp, yp, zp = 0 : periodic
       DEALLOCATE(lxf)
       DEALLOCATE(lxs)
       IF(xp.eq.0) THEN
         DEALLOCATE(wxf)
         DEALLOCATE(wxs)
       ENDIF

       DEALLOCATE(lyf)
       DEALLOCATE(lys)
       IF(yp.eq.0) THEN
         DEALLOCATE(wyf)
         DEALLOCATE(wys)
       ENDIF

       DEALLOCATE(lzf)
       DEALLOCATE(lzs)
       IF(zp.eq.0) THEN
         DEALLOCATE(wzf)
         DEALLOCATE(wzs)
       ENDIF

     END SUBROUTINE ludcmp_deallocate

     SUBROUTINE ludcmp_testalloc

       IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated"
       IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated"
       IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated"
       IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated"

       IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated"
       IF (.not. ALLOCATED(lys)) print *, "lys not allocated"
       IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated"
       IF (.not. ALLOCATED(wys)) print *, "wys not allocated"

       IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated"
       IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated"
       IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated"
       IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated"

     END SUBROUTINE ludcmp_testalloc

     SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp)
      INTEGER, INTENT(IN) :: nx,ny,nz
      INTEGER, INTENT(IN) :: xp,yp,zp
      INTEGER :: ierr

       nxc=nx
       nyc=ny
       nzc=nz

       ! CALL ludcmp_testalloc

!       IF(nyc /= nzc) PRINT*,'ny should be equal nz'

!      xp, yp, zp = 0 : periodic
       IF(xp.eq.0) THEN
         CALL p_lud(1,nxc)
       ELSE
         CALL nonp_lud(1,nxc)
       ENDIF

       IF(yp.eq.0) THEN
         CALL p_lud(2,nyc)
       ELSE
         call nonp_lud(2,nyc)
       ENDIF

       IF(zp.eq.0) THEN
         CALL p_lud(3,nzc)
       ELSE 
         call nonp_lud(3,nzc)
       ENDIF

     END SUBROUTINE ludcmp_calculate


     SUBROUTINE test_nonp_lud1(xx, coef)
      INTEGER :: xx
      REAL(KIND=8), DIMENSION(xx) :: aa
      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
       aa=3.
       aa(1)=0.5 ; aa(2)=4.
       aa(xx-1)=4. ; aa(xx)=0.5

       CALL stdlu(aa,xx,coef)
     END SUBROUTINE test_nonp_lud1


     SUBROUTINE test_nonp_lud2(xx, coef)
      INTEGER :: xx
      REAL(KIND=8), DIMENSION(xx) :: aa
      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef
       aa=5.5
       aa(1)=2./11. ; aa(2)=10.
       aa(xx-1)=10. ; aa(xx)=2./11.

       CALL stdlu(aa,xx,coef)
     END SUBROUTINE test_nonp_lud2


     SUBROUTINE test_p_lud1(xx, coef1, coef2)
      INTEGER :: xx
      REAL(KIND=8) :: a
      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
      a=3. ! first derivative

      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
     END SUBROUTINE test_p_lud1


     SUBROUTINE test_p_lud2(xx, coef1, coef2)
      INTEGER :: xx
      REAL(KIND=8) :: a
      REAL(KIND=8), DIMENSION(xx), INTENT(OUT) :: coef1, coef2
      a=11./2. ! second derivative

      CALL ptdlu(a,xx,coef1,coef2) ! x-direction
     END SUBROUTINE test_p_lud2


     SUBROUTINE nonp_lud(xyz,xx)
      INTEGER :: i,xyz,xx
      REAL(KIND=8), DIMENSION(xx) :: aa
       aa=3.
       aa(1)=0.5 ; aa(2)=4.
       aa(xx-1)=4. ; aa(xx)=0.5
! first derivative
       IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction
       IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction
       IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction
       aa=5.5
       aa(1)=2./11. ; aa(2)=10.
       aa(xx-1)=10. ; aa(xx)=2./11.
! second derivative
       IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction
       IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction
       IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction
     END SUBROUTINE nonp_lud


     SUBROUTINE p_lud(xyz,xx)
      INTEGER :: i,xyz,xx
      REAL(KIND=8) :: a
      a=3. ! first derivative
      IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction
      IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction
      IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction
      a=11./2. ! second derivative
      IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction
      IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction
      IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction
     END SUBROUTINE p_lud

     SUBROUTINE stdlu(a,n,l)
      INTEGER :: n
      REAL(KIND=8), INTENT(IN)  :: a(n)
      REAL(KIND=8), INTENT(OUT) :: l(n)
      REAL(KIND=8) :: d
      INTEGER :: i
      l(1)=1.0d0/a(1)
      DO i=2,n
        d=a(i)-l(i-1)
        l(i)=1.0d0/d
      ENDDO
     END SUBROUTINE stdlu

     SUBROUTINE ptdlu(a,n,l,w)
      INTEGER :: n
      REAL(KIND=8), INTENT(IN) :: a
      REAL(KIND=8), INTENT(OUT) :: l(n),w(n)
      INTEGER :: i
      REAL(KIND=8) :: aa(n),d

      DO i=1,n-1
        aa(i)=a
      ENDDO
      i=n-1
      call stdlu(aa,i,l)
      w(1)=1.0
      DO i=2,n-2
        w(i)=-l(i-1)*w(i-1)
      ENDDO
      w(n-1)=1.0-l(n-2)*w(n-2)
      DO i=1,n-1
        w(i)=w(i)*l(i)
      ENDDO
      d=a
      DO i=1,n-1
        d=d-w(i)*w(i)/l(i)
      ENDDO
      l(n)=1./d
     END SUBROUTINE ptdlu


     SUBROUTINE rhs1np(n,h,x,dx,nd)
      INTEGER,INTENT(IN) :: n,nd
      REAL(KIND=8),INTENT(IN) :: h
      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4

      h1=1.d0/h

      r1=7.d0/3.d0
      r2=1.d0/12.d0
      r3=3.
      a=-1.25
      b=1.
      c=0.25

      DO j=1,nd
        dx(j,n-1)=x(j,n)-x(j,n-2)
        dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2))
        dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3))
        dx(j,2)=x(j,3)-x(j,1)
        IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0.
        IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0.
        dx(j,n-1)=dx(j,n-1)*h1*r3
        dx(j,n)=dx(j,n)*h1
        dx(j,1)=dx(j,1)*h1
        dx(j,2)=dx(j,2)*h1*r3
      ENDDO

      DO i=3,n-2
       DO j=1,nd
        t1=x(j,i+1)-x(j,i-1)
        t2=x(j,i+2)-x(j,i-2)
        dx(j,i)=h1*(r1*t1+r2*t2)
       ENDDO
      ENDDO

     END SUBROUTINE rhs1np


     SUBROUTINE dfnonp(n,h,x,dx,nd,dir)
      INTEGER,INTENT(IN) :: n,nd,dir
      REAL(KIND=8),INTENT(IN) :: h
      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL(KIND=8) :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4

      CALL rhs1np (n,h,x,dx,nd)

      IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction
      IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction
      IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction
     END SUBROUTINE dfnonp

     SUBROUTINE dfp(n,h,x,dx,nd,dir)
      INTEGER,INTENT(IN) :: n,nd,dir
      REAL(KIND=8),INTENT(IN) :: h
      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL(KIND=8) :: r1,r2,h1


      ! print *, "dfnonp received (nd,n)", nd, n

      h1=1./h
      r1=7./3.
      r2=1./12.

      DO j=1,nd
        dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3)))
        dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2)))
        dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1)))
        dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n)))
        dx(j,n-1)=dx(j,n-1)*h1
        dx(j,n)=dx(j,n)*h1
        dx(j,1)=dx(j,1)*h1
        dx(j,2)=dx(j,2)*h1
      ENDDO
 
      DO i=3,n-2
       DO j=1,nd
        dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))
        dx(j,i)=dx(j,i)*h1
       ENDDO
      ENDDO

      IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction
      IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction
      IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction

     END SUBROUTINE dfp

     SUBROUTINE ptdslv(r,n,l,w,nd)
      INTEGER,INTENT(IN) :: n,nd
      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l,w
      INTEGER i,j
      REAL(KIND=8), DIMENSION(nd) :: sum
      DO j=1,nd
        sum(j)=w(1)*r(j,1)
        r(j,1)=r(j,1)*l(1)
      ENDDO
      DO i=2,n-1
       DO j=1,nd
        r(j,i)=r(j,i)-r(j,i-1)
        sum(j)=sum(j)+w(i)*r(j,i)
        r(j,i)=r(j,i)*l(i)
       ENDDO
      ENDDO
      DO j=1,nd
        r(j,n)=l(n)*(r(j,n)-sum(j))
        r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
      ENDDO
      DO i=n-2,1,-1
       DO j=1,nd
        r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
       ENDDO
      ENDDO
     END SUBROUTINE ptdslv

     SUBROUTINE d2fp(n,h,x,dx,nd,dir)
      INTEGER,INTENT(IN) :: n,nd,dir
      REAL(KIND=8),INTENT(IN) :: h
      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL(KIND=8) :: h2,r1,r2,t1,t2


      h2=1./(h*h)
      r1=6.
      r2=3./8.
      DO j=1,nd
        t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2))
        t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3))
        IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
        IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
        dx(j,n-1)=(r1*t1+r2*t2)
        
        t1 = (x(j,1)-2.*x(j,n)+x(j,n-1))
        t2 = (x(j,2)-2.*x(j,n)+x(j,n-2))
        IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
        IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
!        dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
!               +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
        dx(j,n)=(r1*t1+r2*t2)

        t1 = (x(j,2)-2.*x(j,1)+x(j,n))
        t2 = (x(j,3)-2.*x(j,1)+x(j,n-1))
        IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0.
        IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0.       
!        dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
!               +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
        dx(j,1)=(r1*t1+r2*t2)

        t1 = (x(j,3)-2.*x(j,2)+x(j,1))
        t2 = (x(j,4)-2.*x(j,2)+x(j,n))
        IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
        IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0.
!        dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
!               +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
        dx(j,2)=(r1*t1+r2*t2)

        dx(j,n-1)=dx(j,n-1)*h2
        dx(j,n)=dx(j,n)*h2
        dx(j,1)=dx(j,1)*h2
        dx(j,2)=dx(j,2)*h2
      ENDDO
      DO i=3,n-2
       DO j=1,nd
        t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
        t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
        IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
        IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
!        dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
!               +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
        dx(j,i)=(r1*t1+r2*t2)
        dx(j,i)=dx(j,i)*h2
       ENDDO
      ENDDO
      IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction
      IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction
      IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
     END SUBROUTINE d2fp

     SUBROUTINE tdslv(r,n,l,nd)
      INTEGER,INTENT(IN) :: n,nd
      REAL(KIND=8),INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL(KIND=8),INTENT(IN),DIMENSION(n) :: l
      INTEGER i,j
      REAL(KIND=8) t1
      DO j=1,nd
        r(j,1)=r(j,1)*l(1)
      ENDDO
      DO i=2,n
       DO j=1,nd
        t1=r(j,i)-r(j,i-1)
        r(j,i)=l(i)*t1
       ENDDO
      ENDDO
      DO i=n-1,1,-1
       DO j=1,nd
        r(j,i)=r(j,i)-l(i)*r(j,i+1)
       ENDDO
      ENDDO
     END SUBROUTINE tdslv

     SUBROUTINE d2fnonp(n,h,x,dx,nd,dir)
      INTEGER,INTENT(IN) :: n,nd,dir
      REAL(KIND=8),INTENT(IN) :: h
      REAL(KIND=8),INTENT(IN),DIMENSION(nd,n) :: x
      REAL(KIND=8),INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL(KIND=8) :: h2,r1,r2,r3,a,b,c,e,t1,t2


      h2=1./(h*h)
      r1=6.
      r2=3./8.
      r3=12.
      a=13./11.
      b=-27./11.
      c=15./11.
      e=-1./11.

      DO j=1,nd
        dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4))
        dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1))
        dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2))
        dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3))
        IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0.
        IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0.
        IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0.
        IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0.
        dx(j,1)=dx(j,1)*h2
        dx(j,2)=dx(j,2)*h2*r3
        dx(j,n-1)=dx(j,n-1)*h2*r3
        dx(j,n)=dx(j,n)*h2
      ENDDO
      DO i=3,n-2
       DO j=1,nd
        t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
        t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
        IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
        IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.      
 
!        dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
!               +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
        dx(j,i)=(r1*t1+r2*t2)
        dx(j,i)=dx(j,i)*h2
       ENDDO
      ENDDO

      IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction
      IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction
      IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction

     END SUBROUTINE d2fnonp

    END MODULE Compact