subroutine qvhqhv(xh,yh,zh,xv,yv,zv,q11,q12,q22,qvh,qhv, &
            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,4,2,2)::qvh,qhv
       real::mfac,dprodc
       real,dimension(3)::kdir
       real,dimension(0:im+1,0:jm+1,nm,3)::a1h,a2h,a1,a2
       integer::i,j,n,k
       integer,dimension(4)::i1,j1
       integer::isCorner

      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 
      enddo
      enddo
      enddo

      do n=1,nm
      Do j = 1, jm  
      Do i = 1, im  
        k=isCorner(i,j,n,im,jm,nm)
       if(k.eq.0)then
        a1h(i,j,n,1)=(XV(i,j-1,n)-XV(i-1,j,n)) 
        a1h(i,j,n,2)=(YV(i,j-1,n)-YV(i-1,j,n)) 
        a1h(i,j,n,3)=(ZV(i,j-1,n)-ZV(i-1,j,n)) 
       else if(k.ne.2)then
        a1h(i,j,n,1)=(XV(i,j-1,n)-Xh(i,j,n)) 
        a1h(i,j,n,2)=(YV(i,j-1,n)-Yh(i,j,n)) 
        a1h(i,j,n,3)=(ZV(i,j-1,n)-Zh(i,j,n)) 
       else
        a1h(i,j,n,1)=(XV(i-1,j,n)-Xh(i,j,n)) 
        a1h(i,j,n,2)=(YV(i-1,j,n)-Yh(i,j,n)) 
        a1h(i,j,n,3)=(ZV(i-1,j,n)-Zh(i,j,n)) 
       endif
  
       kdir(1)=xh(i,j,n)
       kdir(2)=yh(i,j,n)
       kdir(3)=zh(i,j,n)

       call cprod(kdir,a1h(i,j,n,:),a2h(i,j,n,:))
       call normalize(a1h(i,j,n,:))
       call normalize(a2h(i,j,n,:))
      enddo
      enddo
      enddo
      
      do n=1,nm
      do i=1,im
      do j=1,jm
        i1=(/i-1,i,i-1,i/)
	j1=(/j-1,j-1,j,j/)

	do k=1,4
         qvh(i,j,n,k,1,1)=q11(i1(k),j1(k),n)*dprodc(a1(i1(k),j1(k),n,:),a1h(i,j,n,:)) &
	               +q12(i1(k),j1(k),n)*dprodc(a2(i1(k),j1(k),n,:),a1h(i,j,n,:))
         qvh(i,j,n,k,1,2)=q12(i1(k),j1(k),n)*dprodc(a1(i1(k),j1(k),n,:),a1h(i,j,n,:)) &
	               +q22(i1(k),j1(k),n)*dprodc(a2(i1(k),j1(k),n,:),a1h(i,j,n,:))
         qvh(i,j,n,k,2,1)=q11(i1(k),j1(k),n)*dprodc(a1(i1(k),j1(k),n,:),a2h(i,j,n,:)) &
	               +q12(i1(k),j1(k),n)*dprodc(a2(i1(k),j1(k),n,:),a2h(i,j,n,:))
         qvh(i,j,n,k,2,2)=q12(i1(k),j1(k),n)*dprodc(a1(i1(k),j1(k),n,:),a2h(i,j,n,:)) &
	               +q22(i1(k),j1(k),n)*dprodc(a2(i1(k),j1(k),n,:),a2h(i,j,n,:))
        enddo
       enddo
       enddo
       enddo

      do n=1,nm
      do i=1,im-1
      do j=1,jm-1
        i1=(/i,i+1,i,i+1/)
	j1=(/j,j,j+1,j+1/)

	do k=1,4
         qhv(i,j,n,k,1,1)=dprodc(a1h(i1(k),j1(k),n,:),a1(i,j,n,:))
	 qhv(i,j,n,k,1,2)=dprodc(a2h(i1(k),j1(k),n,:),a1(i,j,n,:))
         qhv(i,j,n,k,2,1)=dprodc(a1h(i1(k),j1(k),n,:),a2(i,j,n,:)) 
         qhv(i,j,n,k,2,2)=dprodc(a2h(i1(k),j1(k),n,:),a2(i,j,n,:))
        enddo
       enddo
       enddo
       enddo
 
      end subroutine qvhqhv