subroutine qintc(xh,yh,zh,xv,yv,zv,q11,q12,q22,im,jm,nm,mfac) 

       implicit none
       integer::im,jm,nm
       real,dimension(0:im+1,0:jm+1,nm)::xh,yh,zh,xv,yv,zv,q11,q12,q22
       real,dimension(0:im+1,0:jm+1,nm,3)::a1h,a2h,a1v1,a1v2,a2v1,a2v2
       real,dimension(0:im+1,0:jm+1,nm,3)::a1,a2,a1n,a2n
       real,dimension(0:im+1,0:jm+1,nm,3)::a1hu,a2hu,a1v1u,a1v2u,a2v1u &
                                      ,a2v2u 
       real::mfac,kdir(3),c1,c2,p11,p12,p21,p22,dprodc
       real,dimension(0:im+1,0:jm+1,nm,8,2,2)::qintc1,qintc2
       integer,dimension(4)::i1,j1
       integer::i,j,n,k,k2
       
      Do n = 1, nm  
      Do j = 0, jm  
      Do i = 0, im  
        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,:))
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 1, jm  
      Do i = 1, im  
        a1h(i,j,n,1) =(Xv(i,j-1,n)-Xv(i-1,j,n))*mfac 
        a1h(i,j,n,2) =(Yv(i,j-1,n)-Yv(i-1,j,n))*mfac 
        a1h(i,j,n,3) =(Zv(i,j-1,n)-Zv(i-1,j,n))*mfac 
        a2h(i,j,n,1) =(Xv(i,j,n)-Xv(i-1,j-1,n))*mfac 
        a2h(i,j,n,2) =(Yv(i,j,n)-Yv(i-1,j-1,n))*mfac 
        a2h(i,j,n,3) =(Zv(i,j,n)-Zv(i-1,j-1,n))*mfac 
        call normalize(a1h(i,j,n,:))
        call normalize(a2h(i,j,n,:))
	kdir(1)=xh(i,j,n)
	kdir(2)=yh(i,j,n)
	kdir(3)=zh(i,j,n)
        call cprod(a2h(i,j,n,:),kdir,a1hu(i,j,n,:))
        call cprod(kdir,a1h(i,j,n,:),a2hu(i,j,n,:))
	c1=dprodc(a1hu(i,j,n,:),a1h(i,j,n,:))
	c2=dprodc(a2hu(i,j,n,:),a2h(i,j,n,:))
	a1hu(i,j,n,:)=a1hu(i,j,n,:)/c1
	a2hu(i,j,n,:)=a2hu(i,j,n,:)/c2
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 0, jm  
      Do i = 0, im
        i1=(/i,i+1,i,i+1/)
	j1=(/j,j,j+1,j+1/)
        do k=1,4
	 p11=dprodc(a1(i,j,n,:),a1h(i1(k),j1(k),n,:)) 
	 p12=dprodc(a1(i,j,n,:),a2h(i1(k),j1(k),n,:)) 
	 p21=dprodc(a2(i,j,n,:),a1h(i1(k),j1(k),n,:)) 
	 p22=dprodc(a2(i,j,n,:),a2h(i1(k),j1(k),n,:)) 
	 qintc1(i,j,n,k,1,1)=q11(i,j,n)*p11+q12(i,j,n)*p21
	 qintc1(i,j,n,k,1,2)=q12(i,j,n)*p11+q22(i,j,n)*p21
	 qintc1(i,j,n,k,2,1)=q11(i,j,n)*p12+q12(i,j,n)*p22
	 qintc1(i,j,n,k,2,2)=q12(i,j,n)*p12+q22(i,j,n)*p22

	 qintc2(i,j,n,k,1,1)=dprodc(a1hu(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k,1,2)=dprodc(a2hu(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k,2,1)=dprodc(a1hu(i1(k),j1(k),n,:),a2(i,j,n,:))
	 qintc2(i,j,n,k,2,2)=dprodc(a2hu(i1(k),j1(k),n,:),a2(i,j,n,:))
        enddo
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 1, jm  
      Do i = 1, im  
        a1v1(i,j,n,:) =a1n(i,j,n,:)+a1n(i-1,j,n,:)
        a2v1(i,j,n,:) =a2n(i,j,n,:)+a2n(i-1,j,n,:) 
        call normalize(a1v1(i,j,n,:))
        call normalize(a2v1(i,j,n,:))
	kdir(1)=0.5*(xv(i,j,n)+xv(i-1,j,n))
	kdir(2)=0.5*(yv(i,j,n)+yv(i-1,j,n))
	kdir(3)=0.5*(zv(i,j,n)+zv(i-1,j,n))
        call cprod(a2v1(i,j,n,:),kdir,a1v1u(i,j,n,:))
        call cprod(kdir,a1v1(i,j,n,:),a2v1u(i,j,n,:))
	c1=dprodc(a1v1u(i,j,n,:),a1v1(i,j,n,:))
	c2=dprodc(a2v1u(i,j,n,:),a2v1(i,j,n,:))
	a1v1u(i,j,n,:)=a1v1u(i,j,n,:)/c1
	a2v1u(i,j,n,:)=a2v1u(i,j,n,:)/c2
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 1, jm  
      Do i = 1, im  
        a1v2(i,j,n,:) =a1n(i,j,n,:)+a1n(i,j-1,n,:)
        a2v2(i,j,n,:) =a2n(i,j,n,:)+a2n(i,j-1,n,:) 
        call normalize(a1v2(i,j,n,:))
        call normalize(a2v2(i,j,n,:))
	kdir(1)=0.5*(xv(i,j,n)+xv(i,j-1,n))
	kdir(2)=0.5*(yv(i,j,n)+yv(i,j-1,n))
	kdir(3)=0.5*(zv(i,j,n)+zv(i,j-1,n))
        call cprod(a2v2(i,j,n,:),kdir,a1v2u(i,j,n,:))
        call cprod(kdir,a1v2(i,j,n,:),a2v2u(i,j,n,:))
	c1=dprodc(a1v2u(i,j,n,:),a1v2(i,j,n,:))
	c2=dprodc(a2v2u(i,j,n,:),a2v2(i,j,n,:))
	a1v2u(i,j,n,:)=a1v2u(i,j,n,:)/c1
	a2v2u(i,j,n,:)=a2v2u(i,j,n,:)/c2
      enddo
      enddo
      enddo

      Do n = 1, nm  
      Do j = 0, jm  
      Do i = 0, im
        i1=(/i,i+1,i,i/)
	j1=(/j,j,j,j+1/)
        do k=1,2
	 k2=k+4
	 p11=dprodc(a1(i,j,n,:),a1v1(i1(k),j1(k),n,:)) 
	 p12=dprodc(a1(i,j,n,:),a2v1(i1(k),j1(k),n,:)) 
	 p21=dprodc(a2(i,j,n,:),a1v1(i1(k),j1(k),n,:)) 
	 p22=dprodc(a2(i,j,n,:),a2v1(i1(k),j1(k),n,:)) 
	 qintc1(i,j,n,k2,1,1)=q11(i,j,n)*p11+q12(i,j,n)*p21
	 qintc1(i,j,n,k2,1,2)=q12(i,j,n)*p11+q22(i,j,n)*p21
	 qintc1(i,j,n,k2,2,1)=q11(i,j,n)*p12+q12(i,j,n)*p22
	 qintc1(i,j,n,k2,2,2)=q12(i,j,n)*p12+q22(i,j,n)*p22

	 qintc2(i,j,n,k2,1,1)=dprodc(a1v1u(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k2,1,2)=dprodc(a2v1u(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k2,2,1)=dprodc(a1v1u(i1(k),j1(k),n,:),a2(i,j,n,:))
	 qintc2(i,j,n,k2,2,2)=dprodc(a2v1u(i1(k),j1(k),n,:),a2(i,j,n,:))
        enddo

        do k=3,4
	 k2=k+4
	 p11=dprodc(a1(i,j,n,:),a1v2(i1(k),j1(k),n,:)) 
	 p12=dprodc(a1(i,j,n,:),a2v2(i1(k),j1(k),n,:)) 
	 p21=dprodc(a2(i,j,n,:),a1v2(i1(k),j1(k),n,:)) 
	 p22=dprodc(a2(i,j,n,:),a2v2(i1(k),j1(k),n,:)) 
	 qintc1(i,j,n,k2,1,1)=q11(i,j,n)*p11+q12(i,j,n)*p21
	 qintc1(i,j,n,k2,1,2)=q12(i,j,n)*p11+q22(i,j,n)*p21
	 qintc1(i,j,n,k2,2,1)=q11(i,j,n)*p12+q12(i,j,n)*p22
	 qintc1(i,j,n,k2,2,2)=q12(i,j,n)*p12+q22(i,j,n)*p22

	 qintc2(i,j,n,k2,1,1)=dprodc(a1v2u(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k2,1,2)=dprodc(a2v2u(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qintc2(i,j,n,k2,2,1)=dprodc(a1v2u(i1(k),j1(k),n,:),a2(i,j,n,:))
	 qintc2(i,j,n,k2,2,2)=dprodc(a2v2u(i1(k),j1(k),n,:),a2(i,j,n,:))
        enddo
      enddo
      enddo
      enddo

      open(21,file='../data_in/grid/qintc.dat',form='unformatted')
      write(21)qintc1
      write(21)qintc2
      close(21)

      end subroutine qintc