subroutine qintc_2(xh1,yh1,zh1,xv1,yv1,zv1,im,jm,nm,mfac) 

       implicit none
       integer:: im,jm,nm
       real,dimension(0:im+1,0:jm+1,nm):: xh1,yh1,zh1,xv1,yv1,zv1
       real,dimension(-1:im+2,-1:jm+2,nm):: xh,yh,zh,xv,yv,zv
       real,dimension(-1:im+2,-1:jm+2,nm,3):: a1,a2,a1n,a2n,a1k,a2k,xi,yj
       real,dimension(-1:im+2,-1:jm+2,nm,3):: b1,b2 
       real:: mfac,kdir(3),c1(3),dprodc,G,g3(3)
       real,dimension(0:im+1,0:jm+1,nm,13,2,2):: qa
       real,dimension(0:im+1,0:jm+1,nm,2,2):: qb
       integer:: i1,j1,i1arr(4),j1arr(4)
       integer:: i,j,n,k,k1,k2,isCorner

       do n=1,nm
       do j=0,jm+1
       do i=0,im+1
       do k=1,13
       do k1=1,2
       do k2=1,2
         qa(i,j,n,k,k1,k2)=0
         qb(i,j,n,k1,k2)=0
       enddo
       enddo
       enddo
       enddo
       enddo
       enddo
       
       xh(0:im+1,0:jm+1,:)=xh1
       yh(0:im+1,0:jm+1,:)=yh1
       zh(0:im+1,0:jm+1,:)=zh1
       xv(0:im+1,0:jm+1,:)=xv1
       yv(0:im+1,0:jm+1,:)=yv1
       zv(0:im+1,0:jm+1,:)=zv1

       call bocoh2l_cb(xh,im,jm,2)
       call bocoh2l_cb(yh,im,jm,2)
       call bocoh2l_cb(zh,im,jm,2)
       call bocov2l_cb(xv,im,jm,2)
       call bocov2l_cb(yv,im,jm,2)
       call bocov2l_cb(zv,im,jm,2)


      Do n = 1, nm  
      Do j = -1, jm+1 
      Do i = -1, im+1  
        a1(i,j,n,1) =(XH(i+1,j,n)-XH(i,j+1,n))*mfac 
        a1(i,j,n,2) =(YH(i+1,j,n)-YH(i,j+1,n))*mfac 
        a1(i,j,n,3) =(ZH(i+1,j,n)-ZH(i,j+1,n))*mfac 
        a2(i,j,n,1) =(XH(i+1,j+1,n)-XH(i,j,n))*mfac 
        a2(i,j,n,2) =(YH(i+1,j+1,n)-YH(i,j,n))*mfac 
        a2(i,j,n,3) =(ZH(i+1,j+1,n)-ZH(i,j,n))*mfac 
        a1n(i,j,n,:)=a1(i,j,n,:)
        a2n(i,j,n,:)=a2(i,j,n,:)
          call normalize(a1n(i,j,n,:))
          call normalize(a2n(i,j,n,:))
        kdir(1)=xv(i,j,n)
        kdir(2)=yv(i,j,n)
        kdir(3)=zv(i,j,n)
          call cprod(a1n(i,j,n,:),kdir,a1k(i,j,n,:))
          call cprod(a2n(i,j,n,:),kdir,a2k(i,j,n,:))
        xi(i,j,n,:)=a1n(i,j,n,:)+a2n(i,j,n,:)+a1k(i,j,n,:)+a2k(i,j,n,:)
        yj(i,j,n,:)=a1n(i,j,n,:)+a2n(i,j,n,:)-a1k(i,j,n,:)-a2k(i,j,n,:)
        c1=a1n(i,j,n,:)+a2n(i,j,n,:)
!################################################################################### 
        if(c1(1)/=0.0.or.c1(2)/=0.0.or.c1(3)/=0.0) then
          xi(i,j,n,:)=xi(i,j,n,:)/(sqrt(dprodc(c1,c1)))*sqrt(2.)*0.5
          yj(i,j,n,:)=yj(i,j,n,:)/(sqrt(dprodc(c1,c1)))*sqrt(2.)*0.5
        else
          print *,'QINTC_2: xi(i,j,n,:),yj(i,j,n,:)=',xi(i,j,n,:),yj(i,j,n,:),i,j,n
        endif
!################################################################################### 

          call cprod(a1(i,j,n,:),a2(i,j,n,:),g3)
        G=dprodc(kdir,g3)
          call cprod(a2(i,j,n,:),kdir,b1(i,j,n,:))
          call cprod(kdir,a1(i,j,n,:),b2(i,j,n,:))
!################################################################################### 
        if(G.ne.0.) then
          b1(i,j,n,:)=b1(i,j,n,:)/G
          b2(i,j,n,:)=b2(i,j,n,:)/G
        else 
          b1(i,j,n,:)=0.
          b2(i,j,n,:)=0.
        end if
!################################################################################### 
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 0, jm
      Do i = 0, im
       
!===========================       
        qb(i,j,n,1,1)=dprodc(b1(i,j,n,:),xi(i,j,n,:))
        qb(i,j,n,1,2)=dprodc(b1(i,j,n,:),yj(i,j,n,:))
        qb(i,j,n,2,1)=dprodc(b2(i,j,n,:),xi(i,j,n,:))!!!!!!!!!!DRAGAN
        qb(i,j,n,2,2)=dprodc(b2(i,j,n,:),yj(i,j,n,:))
       
!============================       
       
        k=1
        do j1=-1,1
        do i1=-1,1
          qa(i,j,n,k,1,1)=dprodc(a1(i,j,n,:),b1(i+i1,j+j1,n,:))
          qa(i,j,n,k,1,2)=dprodc(a1(i,j,n,:),b2(i+i1,j+j1,n,:))
          qa(i,j,n,k,2,1)=dprodc(a2(i,j,n,:),b1(i+i1,j+j1,n,:))
          qa(i,j,n,k,2,2)=dprodc(a2(i,j,n,:),b2(i+i1,j+j1,n,:))

          k=k+1
        enddo
        enddo

          if(isCorner(i+1,j+1,n,im,jm,nm).eq.1.or.   &
             isCorner(i,j+1,n,im,jm,nm).eq.2.or.       &
             isCorner(i+1,j,n,im,jm,nm).eq.3.or.      &
             isCorner(i,j,n,im,jm,nm).eq.4)then
            qa(i,j,n,:,1,1)=1.
            qa(i,j,n,:,1,2)=0.
            qa(i,j,n,:,2,1)=0.
            qa(i,j,n,:,2,2)=1.
          else if(isCorner(i,j,n,im,jm,nm).eq.1)then
            qa(i,j,n,1,1,1)=1. 
            qa(i,j,n,1,1,2)=0.
            qa(i,j,n,1,2,1)=0.
            qa(i,j,n,1,2,2)=1.
          else if(isCorner(i+1,j,n,im,jm,nm).eq.2)then
            qa(i,j,n,3,1,1)=1. 
            qa(i,j,n,3,1,2)=0.
            qa(i,j,n,3,2,1)=0.
            qa(i,j,n,3,2,2)=1.
          else if(isCorner(i,j+1,n,im,jm,nm).eq.3)then
            qa(i,j,n,7,1,1)=1. 
            qa(i,j,n,7,1,2)=0.
            qa(i,j,n,7,2,1)=0.
            qa(i,j,n,7,2,2)=1.
          else if(isCorner(i+1,j+1,n,im,jm,nm).eq.4)then
            qa(i,j,n,9,1,1)=1. 
            qa(i,j,n,9,1,2)=0.
            qa(i,j,n,9,2,1)=0.
            qa(i,j,n,9,2,2)=1.
          endif
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 1, jm-1
      Do i = 1, im-1
        i1arr=(/i-2,i+2,i,i/)
        j1arr=(/j,j,j-2,j+2/)
        do k=1,4
          k2=k+9
          qa(i,j,n,k2,1,1)=dprodc(a1(i,j,n,:),b1(i1arr(k),j1arr(k),n,:))
          qa(i,j,n,k2,1,2)=dprodc(a1(i,j,n,:),b2(i1arr(k),j1arr(k),n,:))
          qa(i,j,n,k2,2,1)=dprodc(a2(i,j,n,:),b1(i1arr(k),j1arr(k),n,:))
          qa(i,j,n,k2,2,2)=dprodc(a2(i,j,n,:),b2(i1arr(k),j1arr(k),n,:))
        enddo
      enddo
      enddo
      enddo

      open(21,file='qintc_2.dat',form='unformatted')
      write(21)qa
      write(21)qb
      close(21)
      print *,im,jm,nm

      end subroutine qintc_2