MODULE Compact
     IMPLICIT NONE

     REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, &
                                            lyf,lys,wyf,wys, &
                                            lzf,lzs,wzf,wzs

     INTEGER :: nxc,nyc,nzc

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

     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

     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) PRINT*, 'work array for lud allocation failed'
       ALLOCATE(lxs(nxc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
       IF(xp.eq.0) THEN
         ALLOCATE(wxf(nxc),STAT=ierr)
         IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
         ALLOCATE(wxs(nxc),STAT=ierr)
         IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
       ENDIF

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

       ALLOCATE(lzf(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
       ALLOCATE(lzs(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
       IF(zp.eq.0) THEN
         ALLOCATE(wzf(nzc),STAT=ierr)
         IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
         ALLOCATE(wzs(nzc),STAT=ierr)
         IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'
       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(OUT) :: 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