MODULE m_compact

  use m_openmpi

     IMPLICIT NONE
     PRIVATE
     REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & 
                                        lyf,lys,wyf,wys, &
                                        lzf,lzs,wzf,wzs
!                                           lyzf,lyzs,wyzf,wyzs

     REAL*8, DIMENSION(:), ALLOCATABLE :: zz1, zz2, vz1, vz2

     INTEGER :: nxc,nyc,nzc
     REAL*8, PARAMETER :: ezero = 1.0e-14
      
     PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp,pdfp,pd2fp

    CONTAINS

     SUBROUTINE ludcmp(nxx,nyy,nzz,xp,yp,zp)
      INTEGER, INTENT(IN) :: nxx,nyy,nzz
      INTEGER, INTENT(IN) :: xp,yp,zp
      INTEGER :: ierr

       nxc=nxx
       nyc=nyy
       nzc=nzz
!       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'
         CALL p_lud(1,nxc)
       ELSE
         CALL nonp_lud(1,nxc)
       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'
         CALL p_lud(2,nyc)
       ELSE
         call nonp_lud(2,nyc)
       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'
         CALL p_lud(3,nzc)
       ELSE 
         call nonp_lud(3,nzc)
       ENDIF

       ALLOCATE(vz1(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'

       ALLOCATE(zz1(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'

       ALLOCATE(vz2(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'

       ALLOCATE(zz2(nzc),STAT=ierr)
       IF(ierr /= 0) PRINT*, 'work array for lud allocation failed'

         call ctdlu (3.0, nzc, vz1, zz1)

         call ctdlu (11.0/2.0, nzc, vz2, zz2)

!       CALL x_lud
!       CALL yz_lud

     END SUBROUTINE ludcmp

     SUBROUTINE nonp_lud(xyz,xx)
      INTEGER :: i,xyz,xx
      REAL*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*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*8 :: a(n),l(n)
      REAL*8 :: d
      INTEGER :: i
      l(1)=1.0/a(1)
      DO i=2,n
        d=a(i)-l(i-1)
        l(i)=1.0/d
      ENDDO
     END SUBROUTINE stdlu

     SUBROUTINE ptdlu(a,n,l,w)
      INTEGER :: n
      REAL*8 :: a,l(n),w(n)
      INTEGER :: i
      REAL*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 ctdlu(lambda,n,v,z)
      INTEGER, INTENT(IN) :: n
      REAL*8, INTENT(IN)  :: lambda
      REAL*8, INTENT(OUT) :: v(n)
      REAL*8, INTENT(OUT) :: z(n)
      INTEGER :: i
      REAL*8 :: sigma, alpha

      v(1)=0.0
      v(2)=lambda
      DO i=3,n
        v(i)=lambda - 1./v(i-1)
      ENDDO

      if (lambda > 0.0) then
        alpha = (-lambda + sqrt(lambda**2 - 4.0))/2.0
      else
        alpha = (-lambda - sqrt(lambda**2 - 4.0))/2.0
      endif

      sigma = (1.0 + alpha**2) / (lambda * (1.0 - alpha**2) * (1.0 - alpha**n))

      DO i=1,n
        z(i) = sigma * (alpha**(i-1) + alpha**(n-i+1))
      ENDDO
     END SUBROUTINE ctdlu


     SUBROUTINE pdfp(x,xu,xl,dx,h,nd,n,nall,dir)

      INTEGER,INTENT(IN) :: nall,n,nd,dir
      REAL*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,2) :: xu, xl
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx

      INTEGER :: i,j
      REAL*8 :: r1,r2,h1

      integer (kind=MPI_INTEGER_KIND) :: idx
      integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
      logical :: recvLow, recvUpp

      integer (kind=MPI_INTEGER_KIND), parameter :: rhstag = 0

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

      !if (myid.eq.0) write(*,*) "parallel dfp commuication"
      if (myid.eq.master) then

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, numprocs-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, numprocs-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, myid+1,     rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      elseif (myid.eq.numprocs-1) then

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0,      rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, 0,      rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      else

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, myid+1, rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      endif

      call MPI_STARTALL (4, requests, mpi_err)

      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)))*h1
       ENDDO
      ENDDO

      recvLow = .true.
      recvUpp = .true.

      do while (recvLow.or.recvUpp)

          call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)

          if (idx.eq.1) then
              DO j=1,nd
                dx(j,1)  =(r1*( x(j,2)-xl(j,2))  +r2*( x(j,3)-xl(j,1)))  *h1
              ENDDO

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

              recvLow = .false.
          endif

          if (idx.eq.2) then
              DO j=1,nd
                dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
              ENDDO

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

              recvUpp = .false.
          endif
      enddo

      IF (dir.eq.1) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
      IF (dir.eq.2) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction
      IF (dir.eq.3) CALL ctdslv(dx,vz1,zz1,nd,n,nall) ! x-direction

     END SUBROUTINE pdfp


     SUBROUTINE pd2fp(x,xu,xl,dx,h,nd,n,nall,dir)
      INTEGER,INTENT(IN) :: nall,n,nd,dir
      REAL*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*8 :: h2,r1,r2,t1,t2

      REAL*8,DIMENSION(nd,2) :: xu, xl

      integer (kind=MPI_INTEGER_KIND) :: idx
      integer (kind=MPI_INTEGER_KIND),dimension(4) :: requests
      logical :: recvLow, recvUpp

      integer (kind=MPI_INTEGER_KIND), parameter :: rhstag = 0

      h2=1./(h*h)
      r1=6.
      r2=3./8.

      !if (myid.eq.0) write(*,*) "parallel d2fp commuication"
      if (myid.eq.master) then

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, numprocs-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, numprocs-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, myid+1,     rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      elseif (myid.eq.numprocs-1) then

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, 0,      rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, 0,      rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      else

        call MPI_SEND_INIT (x(1,n-1), 2*nd, MPI_REAL8, myid+1, rhstag, &
                            MPI_COMM_TASK, requests(1), mpi_err)

        call MPI_SEND_INIT (x(1,1),   2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(2), mpi_err)

        call MPI_RECV_INIT (xl,       2*nd, MPI_REAL8, myid-1, rhstag, &
                            MPI_COMM_TASK, requests(3), mpi_err)

        call MPI_RECV_INIT (xu,       2*nd, MPI_REAL8, myid+1, rhstag, &
                            MPI_COMM_TASK, requests(4), mpi_err)

      endif

      call MPI_STARTALL (4, requests, mpi_err)

      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))
        dx(j,i)=(r1*t1+r2*t2)*h2

       ENDDO
      ENDDO

      recvLow = .true.
      recvUpp = .true.

      do while (recvLow.or.recvUpp)

          call MPI_WAITANY (2, requests(3:), idx, mpi_status, mpi_err)

          if (idx.eq.1) then

              DO j=1,nd

                t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
                t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
                dx(j,1)=(r1*t1+r2*t2)*h2

              ENDDO

              DO j=1,nd

                t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
                t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
                dx(j,2)=(r1*t1+r2*t2)*h2

              ENDDO

              recvLow = .false.
          endif

          if (idx.eq.2) then

              DO j=1,nd

                t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
                t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
                dx(j,n-1)=(r1*t1+r2*t2)*h2

              ENDDO

              DO j=1,nd

                t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
                t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
                dx(j,n)=(r1*t1+r2*t2)*h2

              ENDDO

              recvUpp = .false.
          endif
      enddo

      IF (dir.eq.1) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! x-direction
      IF (dir.eq.2) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! y-direction
      IF (dir.eq.3) CALL ctdslv(dx,vz2,zz2,nd,n,nall) ! z-direction
     END SUBROUTINE pd2fp


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

      h1=1./h

      r1=7./3.
      r2=1./12.
      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
      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*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*8 :: r1,r2,h1

      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 par_dfp(x,xu,xl,dx,h,nd,n,nall,dir)
      INTEGER,INTENT(IN) :: nall,n,nd,dir
      REAL*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*8 :: r1,r2,h1

      REAL*8,DIMENSION(nd,2) :: xu, xl

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

      !if (myid.eq.0) write(*,*) "parallel dfp commuication"
      if (myid.eq.master) then

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, numprocs-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1,     myid+1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      elseif (myid.eq.numprocs-1) then

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, 0,      0, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0,      myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      else

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1, myid+1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      endif

      DO j=1,nd
        dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
        dx(j,n)=  (r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
        dx(j,1)=  (r1*( x(j,2)-xl(j,2))  +r2*( x(j,3)-xl(j,1)))*h1
        dx(j,2)=  (r1*( x(j,3)- x(j,1))  +r2*( x(j,4)-xl(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)))*h1
       ENDDO
      ENDDO


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

     END SUBROUTINE par_dfp

     SUBROUTINE ptdslv(r,n,l,w,nd)
      INTEGER,INTENT(IN) :: n,nd
      REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL*8,INTENT(IN),DIMENSION(:) :: l,w
      INTEGER i,j
      REAL*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 ctdslv(r,va,za,nd,n,nall)
      INTEGER,PARAMETER :: nb = 256
      INTEGER,PARAMETER :: nb2 = 64
      INTEGER,INTENT(IN) :: n,nd,nall
      REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL*8,INTENT(IN),DIMENSION(nall) :: va
      REAL*8,INTENT(IN),DIMENSION(nall) :: za
      REAL*8,DIMENSION(n) :: v,z
      INTEGER i,j
      INTEGER ii,jj

      REAL*8, DIMENSION(nd) :: x1, x0

      INTEGER npart, nbase, nlow, nupp
      INTEGER :: pid, np
      REAL*8, DIMENSION(nb) :: r0, r1

      pid = myid
      np = numprocs

      npart = nall / np
      nbase = pid * npart
      nlow = pid * npart + 1
      nupp = (pid + 1) * npart

      v = va(nlow:nupp)
      z = za(nlow:nupp)

      if (npart.lt.4) then
        ! assertion fail
      endif

      x1 = 0.0
      x0 = 0.0

      do i=1,n
        x0 = x0 + z(i)*r(:,i)
      enddo

      call MPI_REDUCE(x0, x1, nd, MPI_REAL8, MPI_SUM, master, MPI_COMM_TASK, mpi_err)

      ! first process
      if (myid.eq.0) then

        CALL MPI_SSEND(x1, nd, MPI_REAL8, np-1, 100, MPI_COMM_TASK, mpi_err)

        r(:,1) = x1
        r(:,2) = r(:,2) - x1

      ! last process
      elseif (pid.eq.(np-1)) then

        CALL MPI_RECV(x1, nd, MPI_REAL8, master, 100, MPI_COMM_TASK, mpi_status, mpi_err)

        r(:,n) = r(:,n) - x1

      ! intermediate process
      else

      endif

      ! first process
      if (myid.eq.0) then

          DO jj=1,nd,nb

               DO j=jj,jj+nb-1
                r(j,2)=r(j,2)/v(2)
               ENDDO

              DO i=3,n
               DO j=jj,jj+nb-1
                r(j,i)=(r(j,i)-r(j,i-1))/v(i)
               ENDDO
              ENDDO

              r0 = r(jj:jj+nb-1,n)
              CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, jj, MPI_COMM_TASK, mpi_err)

          ENDDO

          DO jj=1,nd,nb

              CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, nd+jj, MPI_COMM_TASK, mpi_status, mpi_err)

              i=n
               DO j=jj,jj+nb-1
                r(j,i) = r(j,i) - r1(j-jj+1)/v(i)
               ENDDO

              DO i=n-1,2,-1
               DO j=jj,jj+nb-1
                r(j,i) = r(j,i) - r(j,i+1)/v(i)
               ENDDO
              ENDDO

          ENDDO

      ! last process
      elseif (pid.eq.(np-1)) then

          DO jj=1,nd,nb

              CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, jj, MPI_COMM_TASK, mpi_status, mpi_err)

              i=1
               DO j=jj,jj+nb-1
                r(j,i)=(r(j,i) - r1(j-jj+1))/v(i)
               ENDDO

              DO i=2,n
               DO j=jj,jj+nb-1
                r(j,i)=(r(j,i) - r(j,i-1))/v(i)
               ENDDO
              ENDDO

          ENDDO

          DO jj=1,nd,nb

              DO i=n-1,1,-1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-r(j,i+1)/v(i)
               ENDDO
              ENDDO

              r0 = r(jj:jj+nb-1,1)
              CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, nd+jj, MPI_COMM_TASK, mpi_err)

          ENDDO

      ! intermediate process
      else

          DO jj=1,nd,nb

              CALL MPI_RECV(r1, nb, MPI_REAL8, pid-1, jj, MPI_COMM_TASK, mpi_status, mpi_err)

              i=1
               DO j=jj,jj+nb-1
                r(j,i)=(r(j,i)-r1(j-jj+1))/v(i)
               ENDDO

              DO i=2,n
               DO j=jj,jj+nb-1
                r(j,i)=(r(j,i)-r(j,i-1))/v(i)
               ENDDO
              ENDDO

              r0 = r(jj:jj+nb-1,n)
              CALL MPI_SEND(r0, nb, MPI_REAL8, pid+1, jj, MPI_COMM_TASK, mpi_err)

          ENDDO

          DO jj=1,nd,nb

              CALL MPI_RECV(r1, nb, MPI_REAL8, pid+1, nd+jj, MPI_COMM_TASK, mpi_status, mpi_err)

              i=n
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-r1(j-jj+1)/v(i)
               ENDDO

              DO i=n-1,1,-1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-r(j,i+1)/v(i)
               ENDDO
              ENDDO

              r0 = r(jj:jj+nb-1,1)
              CALL MPI_SEND(r0, nb, MPI_REAL8, pid-1, nd+jj, MPI_COMM_TASK, mpi_err)

          ENDDO

      endif

     END SUBROUTINE ctdslv


     SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall)
      INTEGER,PARAMETER :: nb = 64
      INTEGER,INTENT(IN) :: n,nd,nall
      REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL*8,INTENT(IN),DIMENSION(nall) :: la,wa
      REAL*8,DIMENSION(n) :: l,w
      INTEGER i,j
      INTEGER ii,jj
      REAL*8, DIMENSION(nd) :: sum

      INTEGER npart, nbase, nlow, nupp
      INTEGER :: pid, np
      REAL*8, DIMENSION(nb) :: r0, r1, sum0
      REAL*8, DIMENSION(nb,2) :: buf

      pid = myid
      np = numprocs

      npart = nall / np
      nbase = pid * npart
      nlow = pid * npart + 1
      nupp = (pid + 1) * npart

      l = la(nlow:nupp)
      w = wa(nlow:nupp)

      if (npart.lt.4) then
        ! assertion fail
      endif

      call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)

      ! first process
      if (myid.eq.0) then

          DO jj=1,nd,nb
               DO j=jj,jj+nb-1
                sum(j)=w(1)*r(j,1)
                r(j,1)=r(j,1)*l(1)
               ENDDO

              DO i=2,n
               DO j=jj,jj+nb-1
                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

              buf(:,1) = r(jj:jj+nb-1,n)
              buf(:,2) = sum(jj:jj+nb-1)
              CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)

          ENDDO

          DO jj=1,nd,nb

              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
              r0 = buf(:,1)
              r1 = buf(:,2)

              i=n
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
               ENDDO
              DO i=n-1,1,-1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
               ENDDO
              ENDDO

          ENDDO

      ! last process
      elseif (pid.eq.(np-1)) then

          DO jj=1,nd,nb

              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
              r0 = buf(:,1)
              sum0 = buf(:,2)

              i=1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-r0(j-jj+1)
                sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
                r(j,i)=r(j,i)*l(i)
               ENDDO

              DO i=2,n-1
               DO j=jj,jj+nb-1
                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=jj,jj+nb-1
                r(j,n)=l(n)*(r(j,n)-sum(j))
               ENDDO

          ENDDO

          DO jj=1,nd,nb

               DO j=jj,jj+nb-1
                r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
               ENDDO

              DO i=n-2,1,-1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
               ENDDO
              ENDDO

              buf(:,1) = r(jj:jj+nb-1,1)
              buf(:,2) = r(jj:jj+nb-1,n)
              CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)

          ENDDO

      ! intermediate process
      else

          DO jj=1,nd,nb

              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
              r0 = buf(:,1)
              sum0 = buf(:,2)

              i=1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-r0(j-jj+1)
                sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
                r(j,i)=r(j,i)*l(i)
               ENDDO

              DO i=2,n
               DO j=jj,jj+nb-1
                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

              buf(:,1) = r(jj:jj+nb-1,n)
              buf(:,2) = sum(jj:jj+nb-1)
              CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)

          ENDDO

          DO jj=1,nd,nb

              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
              r0 = buf(:,1)
              r1 = buf(:,2)

              i=n
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
               ENDDO
              DO i=n-1,1,-1
               DO j=jj,jj+nb-1
                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
               ENDDO
              ENDDO

              buf(:,1) = r(jj:jj+nb-1,1)
              !buf(:,2) = r1
              CALL MPI_SEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)

          ENDDO

      endif

     END SUBROUTINE par_ptdslv

     SUBROUTINE d2fp(n,h,x,dx,nd,dir)
      INTEGER,INTENT(IN) :: n,nd,dir
      REAL*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*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))
        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))
!        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))
!        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))
!        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))
!        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 par_d2fp(x,xu,xl,dx,h,nd,n,nall,dir)
      INTEGER,INTENT(IN) :: nall,n,nd,dir
      REAL*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*8 :: h2,r1,r2,t1,t2

      REAL*8,DIMENSION(nd,2) :: xu, xl

      h2=1./(h*h)
      r1=6.
      r2=3./8.

      !if (myid.eq.0) write(*,*) "parallel d2fp commuication"
      if (myid.eq.master) then

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, numprocs-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1,     myid+1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      elseif (myid.eq.numprocs-1) then

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, 0,      0, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0,      myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      else

              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1, myid+1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
                              MPI_COMM_TASK, mpi_request, mpi_err)

              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
                              MPI_COMM_TASK, mpi_status, mpi_err)

      endif


      DO j=1,nd

        t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
        t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
        dx(j,n-1)=(r1*t1+r2*t2)*h2

        t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
        t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
!        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)*h2

        t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
        t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
!        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)*h2

        t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
        t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
!        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)*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))
!        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)*h2

       ENDDO
      ENDDO

      IF (dir.eq.1) CALL par_ptdslv(dx,lxs,wxs,nd,n,nall) ! x-direction
      IF (dir.eq.2) CALL par_ptdslv(dx,lys,wys,nd,n,nall) ! y-direction
      IF (dir.eq.3) CALL par_ptdslv(dx,lzs,wzs,nd,n,nall) ! z-direction
     END SUBROUTINE par_d2fp

     SUBROUTINE tdslv(r,n,l,nd)
      INTEGER,INTENT(IN) :: n,nd
      REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
      REAL*8,INTENT(IN),DIMENSION(:) :: l
      INTEGER i,j
      REAL*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*8,INTENT(IN) :: h
      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
      INTEGER :: i,j
      REAL*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))
        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))
 
!        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 m_compact