      subroutine werner(iplotall,n,deps,icode,zmin,zmax,nin,xx,z,zz)
c
c "werner deconvolution"
c
c  originally written by W. Anderson
c revised by r.godson,usgs,denver,colorado 4/78
c many revisions since
c
c Last modification: 11 April 91

c Compilation: /Gt1200, link to plot.lib
c
c &parms:
c
c     filnam= input data file name. In single quotes, e.g: 'myfile.pos'
c
c     ichn  = selected post-file channel for processing (1 to 6, default=1).
c             only needed for post files.
c
c     n     = degree of interference polynomial (n.le.5; default n=2).
c             usually n=1 or 2 will suffice.
c           note: n+2 is the highest degree term in system of n+5 eqs.
c             no interference (i.e. n=0) is a special case with 4eqs
c
c     eps   = accuracy or solution tolerance (default 1.0e-2)
c
c     iskip = initial number of points to skip (default 1, for every point
c             see nskip)
c
c     incr  = increment to skip pts for successive passes (default 1,
c             see nskip)
c
c     nskip = final no. pts to skip (default 1) Note: first pass treats every
c             iskip point, 2nd pass every incr point, 3rd pass every 2*incr
c             point, ...last pass treats every nskip point.
c
c     zmin  = min. z to accept as a depth solution (default 0.0)
c
c     zmax  = max. z to accept as a depth solution (default 10.0)
c              note: in plot z-axis ranges from zmin(top) to zmax(botm.)
c              zmin & zmax are given in same units as x-axis (km).
c
c     ipr  = 1 to print in text output file "werner.out" all solutions
c              accepted.
c          = 0 (default) to suppress all solution printing.
c             (the plot will always contain all accepted solutions).
c
c     xscale = horizontal scale in data units (km) per cm.
c
c     vexag = exageration of z (vertical) axis (default=1).
c
c     iplotr= plotter number for plotting system
c             hp7475=5,cga=8,ega=9,vga=10 (default=10)
c
c     name= type of data: post, grid, or form (for formatted string)
c           in single quotes, e.g:  name='post'
c
c     npts= number of input data values (formatted input only)
c
c     deltx= distance between data values (formatted input only)
c
c     iprof= array containing up to 50 row numbers of gridded data
c            (grid imput only)
c
c &
c
c
c  (note that "time" refers always to "distance".)
c     revx = y to reverse axis(i.e. right to left)
c
c
c      common /magdata/ nin,xmin,xmax,amin,amax,zmin,zmax,
c     &xx(512),z(512),zz(512)
      common /work1/zbuf(512),xdata(512)
      common /work2/aa(512),bb(512)
      common /work3/cc(512),pp(512)
      common /work4/ss(512),dummy(512)
c          double precision id(2),rngxy(4),xx(100),yy(100),fil(10)
      dimension xx(nin),z(nin),zz(nin)
      dimension xpsav(512),ypsav(512)
          double precision x,f,a(100),b(10),x1,xnn,xc,xtest,x2,xf,x0,z0
          character*1 sym
          data dval/1.0e38/,size/.08/,dumm/0./,it/0/,dd/0./
c     & icard/9/,iprint/6/,jprint/8/,ldim1/10001/
c          print *,'werner called'
      npts=nin
      do 10 i=1,npts
      xdata(i)=xx(i)
      zbuf(i)=zz(i)
   10 continue
          jsw=0
c          n=2
c          icode=100*n
          eps=1.0e-2
          iskip=1
          incr=1
          nskip=30
          sym='-'
          if(n.eq.0) nn=4
          if(n.gt.0) nn=n+5
          nn1=nn-1
          nn2=nn-2
          tm1=0.
      xstart=xdata(1)
          tm1=xstart
c
c--perform werner deconvolution for current zbuf block--
c
      write(20,202) n
      write(20,201) '   sheet solutions'
      write(20,201) '          x              z            incr'
      kk=0
          deltx=xdata(2)-xdata(1)
101       kskip=iskip
21        kat=1
c          print *,kskip
          xtest=kskip*deltx
22        j=kat
          x1=(kat-1)*deltx
          xnn=(kat-1+kskip*nn1)*deltx
          xc=0.5d0*(x1+xnn)
          i=1
23        x=(j-1)*deltx-xc
          if(zbuf(j).ge.dval) go to 50
          f=zbuf(j)
          a(i)=1.0d0
          a(i+nn)=x
          if(nn.eq.4) go to 400
          x2=x*x
          do 40 l=3,nn2
          a(i+(l-1)*nn)=x2
40        x2=x2*x
400       a(i+nn2*nn)=f
          xf=x*f
          a(i+nn1*nn)=xf
          b(i)=x*xf
41        i=i+1
          j=j+kskip
          if(i.le.nn) go to 23
c          print *,'calling dgelg'
          call dgelg(b,a,nn,1,eps,ks)
c          print *,'returned from dgelg'
c          print *,'ks = ',ks
          if(ks.eq.-1) go to 50
          x0=0.5d0*b(nn)+xc
c          print *,'x0 = ',x0
          if(x0.lt.0.0d0) go to 50
          x=x0-xc
c          print *,'x = ',x,' xtest = ',xtest
          if(dabs(x).gt.xtest) go to 50
          z0=-b(nn1)-x**2
c          print *,'z0 = ',z0
          if(z0.lt.0.0d0) go to 50
          z0=dsqrt(z0)
c          print *,'z0 = ',z0,' zmin = ',zmin,' zmax = ',zmax
          if(z0.lt.zmin.or.z0.gt.zmax) go to 50
c--(x0,z0) accepted as a solution
c          x0=tmin+x0
42        xp=x0
          yp=z0
          do 43 i=2,npts
          if(xp.gt.xx(i)) go to 43
          go to 44
   43     continue
c   44     print*,nin,i,xx(i),xx(i-1),z(i),z(i-1)
c          yp=yp+z(i-1)+(xp-xx(i-1))*(z(i)-z(i-1))/(xx(i)-xx(i-1))
44        yp=yp+z(i-1)
          if(iplotall.eq.1) call vchar(xp, yp,sym,1,800,size,0.,0.,0.)
          write(20,200) xp,yp,kskip
          kk=kk+1
          xpsav(kk)=xp
          ypsav(kk)=yp
          if((kk+ks).eq.512) go to 70
c          print *,icode
50        kat=kat+kskip
          if(kat.le.npts-nn*kskip+1) go to 22
51        kskip=kskip+incr
          if(kskip.le.nskip) go to 21
c
c--end of plot for this zbuf block--
c
          if(jsw.eq.1) go to 70
c
c         compute horizontal derivative and compute solutions
c
c      print *,'splin1 called'
          call splin1(npts,deltx,dumm,zbuf,aa,bb,cc,it,dd,pp,ss)
c      print *,'returned from splin1'
          ksheet=kk
          sym='*'
          jsw=1
          tm1=tm1+deltx
          tmin=tm1
          npts=npts-1
          do 511 i=1,npts
          zbuf(i)=(3.*cc(i)*deltx+2.*bb(i))*deltx+aa(i)
          xdata(i)=xdata(i+1)
511       continue
      write(20,201) '   contact solutions'
      write(20,201) '          x              z            incr'
          go to 101
70        jsw=0
          sym='-'
c          print*
c          print*,ksheet,kk
      write(20,201) '   average sheet solutions'
          call plotav(deps,sym,size,icode,ksheet,0,xpsav,ypsav)
          sym='*'
      write(20,201) '   average contact solutions'
          call plotav(deps,sym,size,icode,kk,ksheet,xpsav,ypsav)
c          if(name.ne.'airborne') go to 98
998        continue
98         return
200       format(2f15.6,i11)
201       format(a)
202       format('   degree',i2,' interference polynomial')
c
c         messages
c
c9991      write(iprint,69991)
c69991     format('0error--n.lt.0.or.n.gt.5')
c          go to 998
c9992      write(iprint,69992)
c69992     format('0error--tm2.le.tm1 .or. time segment cards not',
c     & ' in ascending order')
c          go to 998
c9993      write(iprint,69993)
c69993     format('0error--ichn.lt.1.or.ichn.gt.nchn')
c          go to 998
c69994     format('0warning--time break at t=',e16.8,
c     & '--segment end set to tm2=',e16.8/)
c          t=tm2+deltx
ccc          go to 20
c9995      write(iprint,1600)
c1600      format(' no rows selected for gridded data')
c          go to 998
c9996      write(iprint,1700)
c1700      format(' npts le 0 for user formatted data')
c          go to 998
c9997      write(iprint,1800)
c1800      format(' deltx le 0 for user formatted data')
c          go to 998
c9998      write(iprint,1900)
c1900      format(' end of file while reading user formatted data',/,
c     & ' npts is probably incorrect')
cc          go to 998
          end
      subroutine dgelg(r,a,m,n,eps,ier)
      dimension a(1),r(1)
      double precision r,a,piv,tb,tol,pivi
      if(m)23,23,1
c     search for greatest element in matrix a
1     ier=0
      piv=0.d0
      mm=m*m
      nm=n*m
      do 3 l=1,mm
      tb=dabs(a(l))
      if(tb-piv)3,3,2
2     piv=tb
      i=l
3     continue
      tol=eps*piv
c     a(i) is pivot element.piv contains the absolute value of a(i)
c
c     start elimination loop
      lst=1
      do 17 k=1,m
c     test on singularity
      if(piv)23,23,4
4     if(ier)7,5,7
5     if(piv-tol)6,6,7
6     ier=k-1
7     pivi=1.d0/a(i)
      j=(i-1)/m
      i=i-j*m-k
      j=j+1-k
c     i+k is row index,j+k is column index of pivot element
      do 8 l=k,nm,m
      ll=l+i
      tb=pivi*r(ll)
      r(ll)=r(l)
8     r(l)=tb
c     is elimination terminated
      if(k-m)9,18,18
c     column interchange in matrix a
9     lend=lst+m-k
      if(j)12,12,10
10    ii=j*m
      do 11 l=lst,lend
      tb=a(l)
      ll=l+ii
      a(l)=a(ll)
11    a(ll)=tb
c     row interchange and pivot row reduction in matrix a
12    do 13 l=lst,mm,m
      ll=l+i
      tb=pivi*a(ll)
      a(ll)=a(l)
13    a(l)=tb
c     save column interchangr information
      a(lst)=j
c     element reduction and next pivot search
      piv=0.d0
      lst=lst+1
      j=0
      do 16 ii=lst,lend
      pivi=-a(ii)
      ist=ii+m
      j=j+1
      do 15 l=ist,mm,m
      ll=l-j
      a(l)=a(l)+pivi*a(ll)
      tb=dabs(a(l))
      if(tb-piv)15,15,14
14    piv=tb
      i=l
15    continue
      do 16 l=k,nm,m
      ll=l+j
16    r(ll)=r(ll)+pivi*r(l)
17    lst=lst+m
c     end of elimination loop
c     back substitution and back interchange
18    if(m-1)23,22,19
19    ist=mm+m
      lst=m+1
      do 21 i=2,m
      ii=lst-i
      ist=ist-lst
      l=ist-m
      l=a(l)+.5d0
      do 21 j=ii,nm,m
      tb=r(j)
      ll=j
      do 20 k=ist,mm,m
      ll=ll+1
20    tb=tb-a(k)*r(ll)
      k=j+l
      r(j)=r(k)
21    r(k)=tb
22    return
c     error return
23    ier=-1
      return
      end
      subroutine splin1(m,h,x,y,a,b,c,it,d,p,s)
c--one dimensional cubic spline coefficient determination.
c
c        by  w.l.anderson, u.s. geological survey, denver, colorado
c
c  parms--- m= number of data points .gt. 2
c           h= equal interval option when h.gt.0. (use dummy x here),
c              unequal intervals if h=0. (x required storage)
c           x= indep.var when h=0. (dim .ge. m).
c           y= dependent variable  (dim .ge. m).
c           a,b,c=coeff.arrays (each dim .ge. m)
c                 results are returned in 1st(m-1) elements of a,b,&c.
c                 also used as work arrays during execution.
c           it= type of boundary condition supplied in d array. use
c              it=1 if 1st derivatives given at end points, or
c              it=0 if 2nd derivatives given at end points.
c           d= boundary array (dim 2) at point 1 and m respectively.
c           p,s= work arrays (each dim=m).
c--error return with m=-(abs(m)) if any parm out of range.
c  the resulting cubic spline is of the form:
c     y=y(i)+a(i)*(x-x(i))+b(i)*(x-x(i))**2+c(i)*(x-x(i))**3
c       for i=1,2,...,m-1
c
c
      real*4  x(1),y(1),a(1),b(1),c(1),d(2),p(1),s(1),mul
      if(it.lt.0.or.it.gt.1.or.h.lt.0..or.m.lt.3) go to 999
      n=m-1
      if(it.eq.0) go to 20
c--1st derivative boundaries given
      ne=n-1
      if(h) 999,11,1
c--equal spacing h .gt. 0. and it=1
    1 hh=3.0/h
      do 2 i=1,ne
      b(i)=4.0
      c(i)=1.0
      a(i)=1.0
    2 p(i)=hh*(y(i+2)-y(i))
      p(1)=p(1)-d(1)
      p(ne)=p(ne)-d(2)
c--solution of tridiagonal matrix eq. of order ne
    3 c(1)=c(1)/b(1)
      p(1)=p(1)/b(1)
      do 4 i=2,ne
      mul=1.0/(b(i)-a(i)*c(i-1))
      c(i)=mul*c(i)
    4 p(i)=mul*(p(i)-a(i)*p(i-1))
c--obtain spline coefficients
      a(ne+it)=p(ne)
      i=ne-1
    5 a(i+it)=p(i)-c(i)*a(i+it+1)
      i=i-1
       if(i.ge.1) go to 5
      if(it.eq.0) go to 6
      a(1)=d(1)
      a(m)=d(2)
    6 if(h.eq.0.) go to 14
      hh=1.0/h
      do 7 i=1,n
      mul=hh*(y(i+1)-y(i))
      b(i)=hh*(3.0*mul-(a(i+1)+2.0*a(i)))
    7 c(i)=hh*hh*(-2.0*mul+a(i+1)+a(i))
      return
c--unequal spacing h=0.. and it=1
   11 do 12 i=1,n
   12 s(i+1)=x(i+1)-x(i)
      do 13 i=1,ne
      b(i)=2.0*(s(i+1)+s(i+2))
      c(i)=s(i+1)
      a(i)=s(i+2)
   13 p(i)=3.0*(s(i+1)**2*(y(i+2)-y(i+1))+s(i+2)**2*(y(i+1)-y(i)))/
     $ (s(i+1)*s(i+2))
      p(1)=p(1)-s(3)*d(1)
      p(ne)=p(ne)-s(n)*d(2)
      go to 3
   14 do 15 i=1,n
      hh=1.0/s(i+1)
      mul=(y(i+1)-y(i))*hh**2
      b(i)=3.0*mul-(a(i+1)+2.0*a(i))*hh
   15 c(i)=-2.0*mul*hh+(a(i+1)+a(i))*hh**2
      return
c--2nd derivative boundaries given
   20 ne=n+1
      if(h) 999,31,21
c--equal spacing h .gt. 0 and it=0
   21 hh=3.0/h
      do 22 i=2,n
      b(i)=4.0
      c(i)=1.0
      a(i)=1.0
   22 p(i)=hh*(y(i+1)-y(i-1))
      b(1)=2.0
      b(ne)=2.0
      c(1)=1.0
      c(ne)=1.0
      a(ne)=1.0
      p(1)=hh*(y(2)-y(1))-0.5*h*d(1)
      p(ne)=hh*(y(m)-y(n))+0.5*h*d(2)
      go to 3
c--unequal spacing h=0 and it=0
   31 do 32 i=1,n
   32 s(i+1)=x(i+1)-x(i)
      n1=n-1
      do 33 i=1,n1
      b(i+1)=2.0*(s(i+1)+s(i+2))
      c(i+1)=s(i+1)
      a(i+1)=s(i+2)
   33 p(i+1)=3.0*(s(i+1)**2*(y(i+2)-y(i+1))+s(i+2)**2*(y(i+1)-y(i)))/
     *     (s(i+1)*s(i+2))
      b(1)=2.0
      b(ne)=2.0
      c(1)=1.0
      c(ne)=1.0
      a(ne)=1.0
      p(1)=3.0*(y(2)-y(1))/s(2)-0.5*s(2)*d(1)
      p(ne)=3.0*(y(m)-y(n))/s(m)+0.5*s(m)*d(2)
      go to 3
  999 m=-iabs(m)
      return
      end
c******************************************************************************
      subroutine plotav(eps,sym,size,ipen,msav,ioff,xin,zin)
      dimension xin(1),zin(1),xsav(512),zsav(512)
      character*1 sym
c
      write(20,200) '          x              z              #'
      do 900 i=1,msav
      j=i+ioff
      xsav(i)=xin(j)
      zsav(i)=zin(j)
900   continue
      nsav=msav-ioff
910   x=xsav(1)
      z=zsav(1)
c      print*,x,z
      nsav=nsav-1
      if(nsav.lt.1) go to 999
      do 915 ii=1,nsav
      xsav(ii)=xsav(ii+1)
      zsav(ii)=zsav(ii+1)
915   continue
      ano=1.
      i=0
920   i=i+1
      if(i.gt.nsav) go to 930
      alpha=(x-xsav(i))**2 + (z-zsav(i))**2
c      d=alpha*conjg(alpha)
      d=sqrt(alpha)*100./abs(z)
      if(d.le.eps) then
c      print*,i,xsav(i),zsav(i)
        x=(ano*x+xsav(i))/(ano+1.)
        z=(ano*z+zsav(i))/(ano+1.)
        ano=ano+1.
c        print*,x,z,ano
        nsav=nsav-1
        if(i.gt.nsav) go to 930
        do 925 ii=i,nsav
        xsav(ii)=xsav(ii+1)
        zsav(ii)=zsav(ii+1)
925     continue
        i=i-1
      endif
      go to 920
930   continue
      if(ano.eq.1.) go to 910
      nn=ano
c      sym=char(nn+48)
c      if(nn.gt.9) sym=char(nn+55)
      call vchar(x,z,sym,1,ipen,size,0.,0.,0.)
      write(20,201) x,z,nn
      go to 910
999   return
200   format(a)
201   format(2f15.6,i11)
      end

