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,k2,isCorner
       
       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_oc(xh,im,jm,2)
       call bocoh2l_oc(yh,im,jm,2)
       call bocoh2l_oc(zh,im,jm,2)
       call bocov2l_oc(xv,im,jm,2)
       call bocov2l_oc(yv,im,jm,2)
       call bocov2l_oc(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,:)
	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

	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(i.eq.-1.and.j.eq.1.and.n.eq.5)then
!	  print *,G,b1(i,j,n,:)
!	  print *,a2(i,j,n,:),kdir
!	  stop
        endif
	b1(i,j,n,:)=b1(i,j,n,:)/G
	b2(i,j,n,:)=b2(i,j,n,:)/G
	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,:))
	qb(i,j,n,2,2)=dprodc(b2(i,j,n,:),yj(i,j,n,:))
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 0, jm
      Do i = 0, im
        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,:))
!	  print *,i,j,n,k,qa(i,j,n,k,1,1),qa(i,j,n,k,1,2),qa(i,j,n,k,2,1)  &
!	        ,qa(i,j,n,k,2,2)
	  k=k+1
        enddo
        enddo
!	print *,i,j,n,qb(i,j,n,1,1),qb(i,j,n,1,2),qb(i,j,n,2,1), &
!	         qb(i,j,n,2,2)
	  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,:))
!	  print *,i,j,n,k,a1(i,j,n,:),b1(i1arr(k),j1arr(k),n,:)
!	print *,i,j,n,qa(i,j,n,k2,1,1),qa(i,j,n,k2,1,2),qa(i,j,n,k2,2,1), &
!	         qa(i,j,n,k2,2,2)
        enddo
      enddo
      enddo
      enddo

      open(21,file='../data_in/grid/qintc_2.dat',form='unformatted')
      write(21)qa
      write(21)qb
      close(21)

      end subroutine qintc_2