subroutine qvhqhv2(xh,yh,zh,xv,yv,zv,q11,q12,q22,qh11,qh12,qh22,  &
              qvh2,qhv2,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)::qvh2,qhv2
       real::mfac,dprodc,gvh,sqhsq,gp2,qs11,qs12,qs22
       real,dimension(0:im+1,0:jm+1,nm,3)::a1h,a2h,a1,a2
       real,dimension(0:im+1,0:jm+1,nm)::qh11,qh12,qh22
       integer::i,j,n,k
       integer,dimension(4)::i1,j1,kp
       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

      a1h=0.
      a2h=0.
      qh11=0.
      qh12=0.
      qh22=0.

      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)) *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

        qs11 = a1h (i,j,n,1)**2+a1h(i,j,n,2)**2+a1h(i,j,n,3)**2
        qs22 = a2h(i,j,n,1)**2+a2h(i,j,n,2)**2+a2h(i,j,n,3)**2 
        qs12 = a1h(i, j, n,1) * a2h(i, j, n,1) + a1h(i, j, n,2) &
           * a2h(i, j, n,2) + a1h(i, j, n,3) * a2h(i, j, n,3)
        SQHsq = qs11 * qs22 - qs12 * qs12  
        Qh11 (i, j, n) = qs22 /sqhsq   
        Qh22 (i, j, n) = qs11 /sqhsq  
        Qh12 (i, j, n) = - qs12  /sqhsq
       endif
  
      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
         qvh2(i,j,n,k,1,1)=qh11(i,j,n)*(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,:)))  &
		       +qh12(i,j,n)*(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,:)))  
         qvh2(i,j,n,k,1,2)=qh11(i,j,n)*(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,:)))  &
		       +qh12(i,j,n)*(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,:)))  
         qvh2(i,j,n,k,2,1)=qh12(i,j,n)*(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,:)))  &
		       +qh22(i,j,n)*(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,:)))  
         qvh2(i,j,n,k,2,2)=qh12(i,j,n)*(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,:)))  &
		       +qh22(i,j,n)*(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/)
	kp=(/4,3,2,1/)

	do k=1,4
	  gp2=qvh2(i1(k),j1(k),n,kp(k),1,1)*qvh2(i1(k),j1(k),n,kp(k),2,2)  &
	     -qvh2(i1(k),j1(k),n,kp(k),1,2)*qvh2(i1(k),j1(k),n,kp(k),2,1)  
          if(gp2.ne.0)then
	    qhv2(i,j,n,k,1,1)=qvh2(i1(k),j1(k),n,kp(k),2,2)/gp2
	    qhv2(i,j,n,k,1,2)=-qvh2(i1(k),j1(k),n,kp(k),1,2)/gp2
	    qhv2(i,j,n,k,2,1)=-qvh2(i1(k),j1(k),n,kp(k),2,1)/gp2
	    qhv2(i,j,n,k,2,2)=qvh2(i1(k),j1(k),n,kp(k),1,1)/gp2
          endif
        enddo
       enddo
       enddo
       enddo
 
      end subroutine qvhqhv2