subroutine sqvmod(xa,ya,za,sqv,dxvdx,dyvdx,dzvdx,  &
          dxvdy,dyvdy,dzvdy,q11,q22,q12,k,mfac)

      implicit none
      include 'param_o.h'

      integer,parameter::kstrd=4, lstrd = kstrd, km = (im - 1) * kstrd+1
      real,dimension(1-kstrd:km+kstrd,1-lstrd:km+lstrd,nm)::xa,ya,za  
      real,dimension(0:im+1, 0:jm+1)::sqv
      real,dimension(0:im+1,0:jm+1,nm)::dxvdx,dyvdx,dzvdx,dxvdy,dyvdy,  &
                dzvdy,q11,q12,q22
      integer::k
      real::mfac,distance,d1,d2,f1,f2,qs11,qs12,qs22,d11,d12
      real,dimension(2,3)::dxmdxc
      real,dimension(3)::p1,p2,p3,p4,pv,p14,p23,pp
      real,dimension(2)::xmp,xm1,xm2
      integer::i,j,n,c1,m2,n2,n1,i1,j1

      do n=1,1
      do i=1,k
      do j=1,k
      do c1=1,1
        if(c1.eq.1)then
	  m2=kstrd*(i-1)+1
	  n2=kstrd*(j-1)+1
	  i1=i
	  j1=j
        else if(c1.eq.2)then
	  m2=km-kstrd*i
	  n2=kstrd*(j-1)+1
	  i1=im-i
	  j1=j
        else if(c1.eq.3)then
	  m2=kstrd*(i-1)+1
	  n2=km-kstrd*j
	  i1=i
	  j1=jm-j
        else
	  m2=km-kstrd*i
	  n2=km-kstrd*j
	  i1=im-i
	  j1=jm-j
        endif

	p1(1)=xa(m2,n2,n)
	p1(2)=ya(m2,n2,n)
	p1(3)=za(m2,n2,n)

	p2(1)=xa(m2+4,n2,n)
	p2(2)=ya(m2+4,n2,n)
	p2(3)=za(m2+4,n2,n)


	p3(1)=xa(m2,n2+4,n)
	p3(2)=ya(m2,n2+4,n)
	p3(3)=za(m2,n2+4,n)


	p4(1)=xa(m2+4,n2+4,n)
	p4(2)=ya(m2+4,n2+4,n)
	p4(3)=za(m2+4,n2+4,n)

	pv(1)=xa(m2+2,n2+2,n)
	pv(2)=ya(m2+2,n2+2,n)
	pv(3)=za(m2+2,n2+2,n)


	call middle(p1(1),p1(2),p1(3),p4(1),p4(2),p4(3), &
	     p14(1),p14(2),p14(3))
        call middle(p2(1),p2(2),p2(3),p3(1),p3(2),p3(3), &
	     p23(1),p23(2),p23(3))
         print *,p14
	 print *,pv
	 print *,p23

        d1=distance(p14(1),p14(2),p14(3),pv(1),pv(2),pv(3))
	d2=distance(p23(1),p23(2),p23(3),pv(1),pv(2),pv(3))

	if(d1.gt.d2)then
	  d11=distance(p1(1),p1(2),p1(3),pv(1),pv(2),pv(3))
	  d12=distance(p4(1),p4(2),p4(3),pv(1),pv(2),pv(3))

	  if(d11.gt.d12)then
	    call extend1(p4,pv,pp)
	    dxvdy(i1,j1,n)=(p4(1)-pp(1))*mfac
	    dyvdy(i1,j1,n)=(p4(2)-pp(2))*mfac
	    dzvdy(i1,j1,n)=(p4(3)-pp(3))*mfac

	    call xctoxm(pp,xmp,dxmdxc,n1)
	    call xctoxm(p4,xm1,dxmdxc,n1)
	    call xctoxm(p1,xm2,dxmdxc,n1)
          else
	    call extend1(p1,pv,pp)
	    dxvdy(i1,j1,n)=(pp(1)-p1(1))*mfac
	    dyvdy(i1,j1,n)=(pp(2)-p1(2))*mfac
	    dzvdy(i1,j1,n)=(pp(3)-p1(3))*mfac

	    call xctoxm(pp,xmp,dxmdxc,n1)
	    call xctoxm(p4,xm2,dxmdxc,n1)
	    call xctoxm(p1,xm1,dxmdxc,n1)
          endif

	   f1=sqrt((xmp(1)-xm1(1))**2+(xmp(2)-xm1(2))**2)
	   f2=sqrt((xm1(1)-xm2(1))**2+(xm1(2)-xm2(2))**2)
	   dxvdy(i1,j1,n)=dxvdy(i1,j1,n)*f2/f1
	   dyvdy(i1,j1,n)=dyvdy(i1,j1,n)*f2/f1
	   dzvdy(i1,j1,n)=dzvdy(i1,j1,n)*f2/f1
        else
	  d11=distance(p2(1),p2(2),p2(3),pv(1),pv(2),pv(3))
	  d12=distance(p3(1),p3(2),p3(3),pv(1),pv(2),pv(3))

	  if(d11.gt.d12)then
	    call extend1(p3,pv,pp)
	    dxvdx(i1,j1,n)=(pp(1)-p3(1))*mfac
	    dyvdx(i1,j1,n)=(pp(2)-p3(2))*mfac
	    dzvdx(i1,j1,n)=(pp(3)-p3(3))*mfac

	    call xctoxm(pp,xmp,dxmdxc,n1)
	    call xctoxm(p3,xm1,dxmdxc,n1)
	    call xctoxm(p2,xm2,dxmdxc,n1)

          else
	    call extend1(p2,pv,pp)
	    dxvdx(i1,j1,n)=(p2(1)-pp(1))*mfac
	    dyvdx(i1,j1,n)=(p2(2)-pp(2))*mfac
	    dzvdx(i1,j1,n)=(p2(3)-pp(3))*mfac
	    call xctoxm(pp,xmp,dxmdxc,n1)
	    call xctoxm(p2,xm1,dxmdxc,n1)
	    call xctoxm(p3,xm2,dxmdxc,n1)
          endif

	  f1=sqrt((xmp(1)-xm1(1))**2+(xmp(2)-xm1(2))**2)
	  f2=sqrt((xm1(1)-xm2(1))**2+(xm1(2)-xm2(2))**2)
	  dxvdx(i1,j1,n)=dxvdx(i1,j1,n)*f2/f1
	  dyvdx(i1,j1,n)=dyvdx(i1,j1,n)*f2/f1
	  dzvdx(i1,j1,n)=dzvdx(i1,j1,n)*f2/f1
        endif

	qs11=dxvdx(i1,j1,n)**2+dyvdx(i1,j1,n)**2+dzvdx(i1,j1,n)**2
	qs22=dxvdy(i1,j1,n)**2+dyvdy(i1,j1,n)**2+dzvdy(i1,j1,n)**2
	qs12=dxvdx(i1,j1,n)*dxvdy(i1,j1,n)+dyvdx(i1,j1,n)*dyvdy(i1,j1,n) &
	     +dzvdx(i1,j1,n)*dzvdy(i1,j1,n)

        if(n.eq.1)then
	  sqv(i1,j1)=sqrt(qs11*qs22-qs12*qs12)
        endif
	q11(i1,j1,n)= qs22/(sqv(i1,j1)**2)
	q22(i1,j1,n)= qs11/(sqv(i1,j1)**2)
	q12(i1,j1,n)=-qs12/(sqv(i1,j1)**2)

	enddo
	enddo
	enddo
	enddo

	end subroutine sqvmod