subroutine dmxinter()

      implicit double precision (a-h,o-z)

#include "divcon.dim"
#include "divcon.h"

      real*8 xiia(3),yiia(3)
      integer n

       n = 3 ! interpolation order
       xiia(1)=1*1.0d0
       xiia(2)=2*1.0d0
       xiia(3)=3*1.0d0
       x=4.0d0

       ! interpolation of the diagonal block
       do i=1,iimat(natoms+1)-1
         yiia(1)=pdiag3(i)
         yiia(2)=pdiag2(i)
         yiia(3)=pdiag1(i)
         call polint(xiia,yiia,n,x,pdiag(i))
       enddo

       ! interpolation of the non diagonal part
       do iatm = 2, natoms
         do ijpair = ip1(iatm), ip1(iatm+1)-1
           jatm = ipair(ijpair)
           ! ijpair = interaction between atom iatm and atom jatm
           ! -> does it exists in previous pairs?

           ! in ipair3
           ijpair3 = -1
           do i = ip13(iatm), ip13(iatm+1)-1
              if (ipair3(i).eq.jatm) then
                ijpair3 = i
                exit
              endif
           end do

           ! in ipair2
           ijpair2 = -1
           do i = ip12(iatm), ip12(iatm+1)-1
              if (ipair2(i).eq.jatm) then
                ijpair2 = i
                exit
              endif
           end do

           ! in ipair1
           ijpair1 = -1
           do i = ip11(iatm), ip11(iatm+1)-1
              if (ipair1(i).eq.jatm) then
                ijpair1 = i
                exit
              endif
           end do

           do ii = 0,ijmat(ijpair+1)-ijmat(ijpair)-1
             if (ijpair3.gt.0) then
               yiia(1)=pdiat3(ijmat3(ijpair3)+ii)
             else
               yiia(1) = 0.0d0
             endif
             if (ijpair2.gt.0) then
               yiia(2)=pdiat2(ijmat2(ijpair2)+ii)
             else
               yiia(2) = 0.0d0
             endif
             if (ijpair1.gt.0) then
               yiia(3)=pdiat1(ijmat1(ijpair1)+ii)
             else
               yiia(3) = 0.0d0
             endif
             call polint(xiia,yiia,n,x,pdiat(ijmat(ijpair)+ii))
           end do
         end do
       end do
       END

c_______________________________________________________________________

      subroutine polint(xa,ya,n,x,y)

      implicit none
      integer i,j,n
      real*8 xa(n),ya(n),x,y,yi
     

       y=0.0d0
      do i=1,n
        yi=ya(i)
        do j=1,n 
          if (j.ne.i) then
          yi=yi*(x-xa(j))/(xa(i)-xa(j)) 
          endif  
        enddo
       y=y+yi
      enddo

      END