subroutine solvlr(a,x,n,lsolv,ks)
*
*     lr-faktorisation,  l*r=a; solves the system lrx=b
*
*     a   full matrix in, lr-faktorised matrix out
*     x   rhs in, solution out
*     n   number of equations
*     lsolv=1 lr-faktorisation only 
*     lsolv=2 backward substituion only (a lr-faktorised)
*     lsolv=3 both lr-faktorisation and back substitution
*     ks=0  lr-faktorisation ok
*     ks=1  a-matrix singular
*
*     nb:
*     l     l-matrix is stored in left triangle of a with supressed
*           1:s in the main diagonal
*     r     r-matrix is stored in right triangle of a
*
      implicit double precision (a-h,o-z)
      dimension a(n,*),x(*)
*
      go to (10,50,10),lsolv
   10 continue
*
*     tolerance of the pivot-element
*
      tol=1.d-20
      ks=0
*
*     for all rows in a
*
         do 40 i=2,n
*
         ii=i-1
*
*        control of the pivot-element
*
         pivot=a(ii,ii)
         if (abs(pivot).gt.tol) go to 20
         ks=1
         go to 100
   20    continue
*
*        blank all elements below the pivot-element
*
            do 40 j=i,n
*
*           compute multiplicator
*
            biga=a(j,ii)/pivot
*
*           subtract row ii*biga from row i
*
               do 30 k=ii,n
               a(j,k)=a(j,k)-a(ii,k)*biga
   30          continue
*
*         save the multiplicator biga
*
            a(j,ii)=biga
   40       continue
*
      if (lsolv.eq.1) go to 100
   50 continue
*
*     forward substitution ly=b
*
         do 70 i=2,n
*
         sum=x(i)
         ii=i-1
            do 60 j=1,ii
            sum=sum-a(i,j)*x(j)
   60       continue
         x(i)=sum
   70    continue
*
*    backward substitution rx=y
*
      x(n)=x(n)/a(n,n)
      kk=n-1
         do 90 k=1,kk
         l=n-k
         sum=x(l)
         ll=l+1
            do 80 m=ll,n
            sum=sum-a(l,m)*x(m)
   80       continue
         x(l)=sum/a(l,l)
   90    continue
*
  100 continue
*
      return
      end
      subroutine backgr(nv,coeg,lquad,rq)
      implicit double precision(a-h,o-z)
*
      dimension coeg(2,*),lquad(7,*),rq(3,*)
*
      xma=-1.d9
      yma=-1.d9
      ymi=1.d9
      xmi=1.d9
      do 100 ip=1,nv
      xt=coeg(1,ip)
      yt=coeg(2,ip)
      xma=dmax1(xma,xt)
      yma=dmax1(yma,yt)
      ymi=dmin1(ymi,yt)
100   xmi=dmin1(xmi,xt)
      rq(1,1)=(xmi+xma)/2.d0
      rq(2,1)=(ymi+yma)/2.d0
      rq(3,1)=dmax1(xma-xmi+0.2d0,yma-ymi+0.2d0)
*
*      nquad=1
*      istart=1
*      call zero(lquad,nv*7*4)
*      do 200 ip=1,nv
*      call addq(lquad,nquad,rq,ip,coeg)
*200   continue
*
      return
      end
      subroutine phsmo(lpoin,lfapo,lface,coeg,ibou,npoin)
      implicit double precision (a-h,o-z)
      dimension lpoin(*),lfapo(3,*),lface(3,*),
     x coeg(2,*),ibou(*)
*
      eps=1.d-9
      do 20 illa=1,10
      do 10 i=1,npoin
      if(ibou(i).ge.0) goto 10
      istart=lpoin(i)
      xm=0.d0
      ym=0.d0
      rnr=0.d0
      xp=coeg(1,i)
      yp=coeg(2,i)
*
90    continue
      itr=lfapo(3,istart)
      i2=itr
      if(itr.lt.0) itr=2
      do 100 ik=1,itr
      iface=lfapo(ik,istart)
      inod=lface(1,iface)
      if(inod.eq.i) inod=lface(2,iface)
      xp=coeg(1,inod)
      yp=coeg(2,inod)
      rnr=rnr+1.d0
      xm=xm+xp
      ym=ym+yp
100   continue
      istart=-i2
      if(i2.lt.0) goto 90
*
*     change coeg, laplacian smoothing
*
      xp=xm/rnr
      yp=ym/rnr
      coeg(1,i)=xp
      coeg(2,i)=yp
10    continue
20    continue
      return
      end




      subroutine checkdir(
     i npoin,coeg,idpnt,
     o area)
      implicit double precision (a-h,o-z)
      dimension coeg(2,*),idpnt(*)
*
*     check direction of the polygon enclosed
*     by the segments joining the points.
*     the area is negative if the segments form 
*     a clockwise polygon
*
      area=0.d0
      do 10 i=1,npoin-1
         ip=idpnt(i)
         ip2=idpnt(i+1)
         area=area+coeg(1,ip)*coeg(2,ip2)-coeg(1,ip2)*coeg(2,ip)
10    continue
      np=idpnt(npoin)
      i1=idpnt(1)
      area=(area+coeg(1,np)*coeg(2,i1)-coeg(1,i1)*coeg(2,np))/2.d0
      return
      end



      subroutine shepard(
     i x,y,npoin,xpoin,ypoin,fpoin,
     o fxy)
      implicit double precision (a-h,o-z)
      dimension xpoin(*),ypoin(*),fpoin(*)
*
*     use shepard's distance-weighted metod
*     for interpolating scattered data.
*     the subroutine returns the interpolated 
*     function fxy evaluated at (x,y). 
*
*     xpoin - x-coordiantes for all points
*     ypoin - y-coordiantes for all points
*     fpoin - function values for all points
*     npoin - number of points
*
      fxy=0.d0
      denom=0.d0
      do 10 i=1,npoin
         xi=xpoin(i)
         yi=ypoin(i)
         dist=dsqrt((x-xi)*(x-xi)+(y-yi)*(y-yi))
         dist=dmax1(dist,1.d-10)
         fxy=fxy+fpoin(i)/dist
         denom=denom+1.d0/dist
10    continue
      fxy=fxy/denom
      return
      end

      subroutine shepard2(
     i x,y,npoin,xpoin,ypoin,fpoin,
     o fxy)
      implicit double precision (a-h,o-z)
      dimension xpoin(*),ypoin(*),fpoin(*)
      dimension e(6,6),v(6)
      data weight/1.d0/
*
*     use shepard's distance-weighted metod
*     for interpolating scattered data.
*     the subroutine returns the interpolated 
*     function fxy evaluated at (x,y). 
*
*     xpoin - x-coordiantes for all points
*     ypoin - y-coordiantes for all points
*     fpoin - function values for all points
*     npoin - number of points
*
      fxy=0.d0
      denom=0.d0
      do 10 i=1,6
         v(i)=0.d0
         do 10 j=1,6
            e(i,j)=0.d0
10    continue
      do 20 i=1,npoin
         xi=xpoin(i)
         yi=ypoin(i)
         x2=xi*xi
         y2=yi*yi
         term=weight*((xi-x)*(xi-x)+(yi-y)*(yi-y))
*
         xterm=xi*term
         yterm=yi*term
         xxterm=x2*term
         yyterm=y2*term
         xyterm=xi*yterm
*
         e(1,1)=e(1,1)+term
         e(1,2)=e(1,2)+xterm
         e(1,3)=e(1,3)+yterm
         e(1,4)=e(1,4)+xyterm
         e(1,5)=e(1,5)+xxterm
         e(1,6)=e(1,6)+yyterm
         e(2,4)=e(2,4)+x2*yterm
         e(2,5)=e(2,5)+x2*xterm
         e(2,6)=e(2,6)+y2*xterm
         e(3,6)=e(3,6)+y2*yterm
         e(4,4)=e(4,4)+x2*yyterm
         e(4,5)=e(4,5)+x2*xyterm
         e(4,6)=e(4,6)+y2*xyterm
         e(5,5)=e(5,5)+x2*xxterm
         e(6,6)=e(6,6)+y2*yyterm
*
         zterm=fpoin(i)*term
         v(1)=v(1)+zterm
         v(2)=v(2)+xi*zterm
         v(3)=v(3)+yi*zterm
         v(4)=v(4)+xi*yi*zterm
         v(5)=v(5)+x2*zterm
         v(6)=v(6)+y2*zterm
20    continue
      e(2,2)=e(1,5)
      e(2,3)=e(1,4)
      e(3,3)=e(1,6)
      e(3,4)=e(2,6)
      e(3,5)=e(2,4)
      e(5,6)=e(4,4)
      do 40 i=1,5
      do 30 j=i+1,6
         e(j,i)=e(i,j)
30    continue
40    continue
      call solvlr(e,v,6,3,ks)
      fxy=v(1)+x*(v(2)+y*v(4)+x*v(5))+y*(v(3)+y*v(6))
      return
      end