c WERNER.FOR
c
c
c  *** uses 'werner deconvolution' on airborne profiles where
c    input may be a post file, a standard grid file or a file containing
c    user formatted data
c    output--plot-- depth solution plots--the number of
c                      plots depends on time selected segments,etc.
c
c--werner.for-- magnetic (profile) interpretation/version 2 --  11/02/73.
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
          common/mintc1/idplot,ind,revx,sym,crt
          common/mintco/zbuf,xdata,npts,iplotr,deltx,
     & xscale,vexag,delx,dely
          character*40 ind,ind1,ind2
          character*40 idplot
          character*80 ifmt
          character*50 fname,filnam,outfile
          character*4 name
        character*8 idpos
        character*1 ff(2)
        integer*2 iff
          dimension zbuf(3000),pos(6)
C          double precision xx(100),yy(100)
          dimension iprof(50),title(14),pgm(2),xdata(3000),
     & aa(3000),bb(3000),cc(3000),pp(3000),ss(3000)
          character*1 sym,star,revx,crt(27000)
c       change namelist to common
          common/parms/ichn,n,eps,iskip,incr,nskip,zmin,zmax,
     & ipr,filnam,ifmt,name,iprof,outfile
          double precision x,f,a(100),b(10),x1,xnn,xc,xtest,x2,xf,x0,z0
        equivalence (iff,ff)
        data iff/#0C00/
          data dval/1.0e38/,size/.08/,dumm/0./,it/0/,dd/0./,
     & icard/9/,iprint/6/,jprint/8/
          in=10                   
          write(iprint,1000)
1000      format(' enter command segment name: ',$)
          read(5,500) fname
500   format(a)
          open(icard,file=fname,status='old')
2         ichn=0
          jsw=0
          n=2
          eps=1.0e-2
          iskip=1
          incr=1
          nskip=30
          zmin=0.0
          zmax=10.0
          xscale=1
          vexag=1.0
          npts=0
          iprof(1)=0
          deltx=0
          ipr=1
          iplotr=10
          sym='x'
          idplot=' '
          filnam=' '
          name=' '
          outfile='werner.out'
          ifmt=' '
          revx='n'
          nsvl=10
          nshl=6
c
c--read &parms
c
        call namemc(icard)
c          read(icard,parms,end=99)
c convert xscale to km/inch, as per original program
          xscale=xscale*2.54
      print*,'Enter plot title:'
      read(5,500)idplot
          print*,'Text output in file "werner.out"'
          open(jprint,file=outfile,status='unknown',form='formatted')
          if(name.eq.'GRID') name='grid'
          if(name.eq.'POST') name='post'
          if(xscale.le.0.) xscale=6.35
          if(vexag.le.0.) vexag=1.
          if(n.lt.0.or.n.gt.5) go to 9991
          if(n.eq.0) nn=4
          if(n.gt.0) nn=n+5
          nn1=nn-1
          nn2=nn-2
          if(iplotr-9) 66,67,68
66        idim=4
          go to 69
67        idim=16
          go to 69
68        idim=18
69        continue
          if(name.eq.'grid'.or.name.eq.'post')
     & open(in,file=filnam,status='old',form='unformatted')
c  select type of input data
          if(name.eq.'post') go to 11
          if(name.eq.'grid') go to 3
      if(name.ne.'form')then
        print*,'Check parameter "name".'
        stop
      endif
c
c  read formatted data
c
          if(npts.le.0) go to 9996
          if(deltx.le.0) go to 9997
          open(in,file=filnam,status='old')
          if(ifmt.eq.' ') then
            read(in,*,end=9998) (zbuf(i),i=1,npts)
          else
            read(in,ifmt,end=9998) (zbuf(i),i=1,npts)
          endif
          xdata(1)=0.
          do 111 i=2,npts
          xdata(i)=xdata(i-1)+deltx
111       continue
          tm1=0.
          tmin=0.
          revx='n'
c
c         plot user formatted input data
c
          write(jprint,1050) ff(2),filnam
1050      format(a1,11x,'input segment - ',a50,/)
          ind1='input data'
          ind2='depths'
112       ind=ind1
          call mintplot(1,zmin,zmax,ksw,idim)
          ind=ind2
          call mintplot(2,zmin,zmax,ksw,idim)
          go to 101
c
c         post-file data
c
11    continue
      i=0
      sum=0
700   read(in,end=701)idpos,xp,yp,pos
      i=i+1
      xdata(i)=sqrt(xp*xp+yp*yp)
      zbuf(i)=pos(ichn)
      if(i.gt.1)sum=sum+xdata(i)-xdata(i-1)
      go to 700
701   npts=i
      if(npts.eq.0)then
        print*,'Abort: problem with post input.'
        stop
      endif
      deltx=sum/npts
      print*,'Average data spacing = ',deltx
      xstart=xdata(1)
c
          write(iprint,704) filnam,title,npts,xstart,deltx
          write(jprint,703) ff(2),filnam,title,npts,xstart,deltx
703     format(a1,1x,'input segment: ',a50,/,
     & ' title = ',14a4,/,' no. of data points = ',i3,/,
     & ' first point coordinate = ',e15.7,/,
     & ' spacing of data points = ', e15.7)
704      format('1input segment: ',a50,/,
     & ' title = ',14a4,/,' no. of data points = ',i3,/,
     & ' first point coordinate = ',e15.7,/,
     & ' spacing of data points = ',e15.7)

          tm1=xstart
          tmin=tm1
c
c         plot input post file
c
          ind1='input post file'
          ind2='depths'

997       ind=ind1
          call mintplot(1,zmin,zmax,ksw,idim)
          ind=ind2
          call mintplot(2,zmin,zmax,ksw,idim)
          go to 101
c  end read post-file data
c  read gridded data
3         ieof=0
          isw=0
          k=0
8         il=1
          k=k+1
          if(isw.ne.0) go to 92
          if(iprof(1).eq.0) go to 9995
          read(in,end=99) title,pgm,ncol,nrow,nz,xstart,deltx,yo,dy
          write(iprint,1450) filnam,title,pgm,nrow,ncol,
     & xstart,deltx,yo,dy
           write(jprint,1451) ff(2),filnam,title,
     & pgm,nrow,ncol,xstart,deltx,yo,dy
1451    format(a1,1x,'input segment: ',a50,/,
     & ' title = ',14a4,/,' creation program id = ',2a4,/,
     & ' no. of rows = ',i3,3x,' no. of columns = ',i3,/,
     & ' first column coordinate = ',e15.7,/,
     & ' spacing of columns = ', e15.7,/,
     & ' first row coordinate = ',e15.7,/,
     & ' spacing of rows = ',e15.7,/)
1450      format('1input segment: ',a50,/,
     & ' title = ',14a4,/,' creation program id = ',2a4,/,
     & ' no. of rows = ',i3,3x,' no. of columns = ',i3,/,
     & ' first column coordinate = ',e15.7,/,
     & ' spacing of columns = ',e15.7,/,
     & ' first row coordinate = ',e15.7,/,
     & ' spacing of rows = ',e15.7,/)
          do 90 i=1,50
          if(iprof(i).eq.0) go to 91
90        continue
          i=51
91        nprof=i-1
          lk=0
92        read(in,end=99) dlt1,(zbuf(j),j=1,ncol)
          lk=lk+1
          if(lk.eq.iprof(k)) go to 93
          go to 92
93        j=0
          tm1=xstart
          do 994 i=1,ncol
          if(zbuf(i).ge.dval) go to 993
          j=j+1
          zbuf(j)=zbuf(i)
          xdata(j)=xstart+(i-1)*deltx
          go to 994
993       if(j.gt.0) go to 994
          tm1=xstart+(i*deltx)
994       continue
          tmin=tm1
          npts=j
          write(iprint,995) iprof(k)
          write(jprint,995) iprof(k)
995       format(/,' processing row ',i5)
c
c         plot input row of gridded data
c
          ind1='input grid row'
          ind2='depths'

          ind=ind1
          call mintplot(1,zmin,zmax,ksw,idim)
          ind=ind2
          call mintplot(2,zmin,zmax,ksw,idim)
          go to 101
c  end read grid data
c
c
c--perform werner deconvolution for current zbuf block--
c
101       kskip=iskip
21        kat=1
          xtest=kskip*deltx
          if(kskip.eq.iskip) go to 211
          write(jprint,621) kskip,xtest
621       format(56x,'iskip=',i2,' xtest=',f9.3)
          go to 22
211       write(jprint,62) ff(2),kskip,xtest
62        format(a1,10x,'x',12x,'xc',10x,'dist',9x,
     & 'depth',3x,'iskip=',i2,
     & ' xtest=',f9.3//)
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
          call dgelg(b,a,nn,1,eps,ks)
          if(ks.eq.-1) go to 50
          x0=0.5d0*b(nn)+xc
          if(x0.lt.0.0d0) go to 50
          x=x0-xc
          if(dabs(x).gt.xtest) go to 50
          z0=-b(nn1)-x**2
          if(z0.lt.0.0d0) go to 50
          z0=dsqrt(z0)
          if(z0.lt.zmin.or.z0.gt.zmax) go to 50
c--(x0,z0) accepted as a solution
          x0=tmin+x0
          if(ipr.eq.0) go to 42
          xc=tmin+xc
          x1=tmin+x1
          write(jprint,64) x1,xc,x0,z0
64        format(1x,4(2x,f11.3))
42        xp=x0
          yp=z0
c
c         plot individual solution positions
c
c          if(ksw.eq.1) call pplot(3,nscale,nhl,nshl,nvl,nsvl,
c     & crt,xmaxp,xminp,ymaxp,yminp,sym,yp,xp,1,1,' ',0)
c          if(ksw.eq.0) call pplot(3,nscale,nhl,nshl,nvl,nsvl,
c     & crt,xmaxp,xminp,ymaxp,yminp,sym,xp,yp,1,1,' ',0)
          call vchar(xp, yp,sym,1,0,size,0.,0.,0.)
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
          call mintplot(3,zmin,zmax,ksw,idim)
          call splin1(npts,deltx,dumm,zbuf,aa,bb,cc,it,dd,pp,ss)
          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
2000      format(a1,'following plots and data result from horizontal',
     & ' derivative data')
          ind1='horizontal derivative'
          ind2='derivative depths'
          write(jprint,2000)ff(2)
          if(name.eq.'grid') go to 997
          go to 112
70        jsw=0
          if(name.eq.'grid') go to 43
          if(name.ne.'airborne') go to 98
          if(ieof.eq.1) idim=2
          call mintplot(3,zmin,zmax,ksw,idim)
          if(ieof.eq.1) go to 99
          if(t.lt.tm2) go to 71
c
c--next segment,if any of airborne data?
c
          revx='n'
          read(icard,*,end=99) t1,t2,irevx
          if(irevx.eq.0) revx='y'
          write(iprint,600) t1,t2,revx
600   format('1',11x,'Selected segment tm1,tm2',2e16.8,5x,
     &'revx =',a1,/)
          if(t2.le.t1.or.t1.lt.tm2) go to 9992
          tm1=t1
          tm2=t2
          tmin=tm1
          isw=0
          go to 8
71        tmin=tmax
          tm1=tmin
          isw=1
          go to 8
c
c         next row of gridded data,if any
c
43        if(k.eq.nprof) idim=2
          call mintplot(3,zmin,zmax,ksw,idim)
          if(k.eq.nprof) go to 99
          isw=1
          go to 8
98        idim=2
          call mintplot(3,zmin,zmax,ksw,idim)
99        write(iprint,1100) outfile
1100      format('0output is in file ',a50)
998       close(jprint)
          close(icard)
          close(in)
          stop
c
c         messages
c
9991      write(iprint,69991)
69991     format('0error--n.lt.0.or.n.gt.5')
          go to 998
9992      write(iprint,69992)
69992     format('0error--tm2.le.tm1 .or. time segment cards not',
     & ' in ascending order')
          go to 998
9993      write(iprint,69993)
69993     format('0error--ichn.lt.1.or.ichn.gt.nchn')
          go to 998
69994     format('0warning--time break at t=',e16.8,
     & '--segment end set to tm2=',e16.8/)
          t=tm2+deltx
cc          go to 20
9995      write(iprint,1600)
1600      format(' no rows selected for gridded data')
          go to 998
9996      write(iprint,1700)
1700      format(' npts le 0 for user formatted data')
          go to 998
9997      write(iprint,1800)
1800      format(' deltx le 0 for user formatted data')
          go to 998
9998      write(iprint,1900)
1900      format(' end of file while reading user formatted data',/,
     & ' npts is probably incorrect')
          go to 998
          end
          subroutine mintplot(ic,ymin,ymax,ksw,idim)
          common/mintc1/idplot,ind,revx,ysym,crt
          common/mintco/zbuf,xdata,npts,iplotr,deltx,
     & xscale,vexag,delx,dely
          character*1 revx,sym,crt(27000),xsym,ysym
          character*40 ind
          character*40 idplot
          dimension xd(2),yd(2),xs(4),ys(4),jdim(18),
     & xdata(3000),zbuf(3000),iddplot(10),indd(3)
          equivalence(iddplot(1),idplot),(indd(1),ind)
        equivalence(sym,iddplot),(xsym,indd)
          data jdim(1)/0/,iprint/6/
c
c         ic = 1 = input profile
c              2 = depth profile
c              3 = close plot system
c
          zmin=ymin
          zmax=ymax
          go to (10,20,70),ic
10        call minmax(zbuf,npts,zmin,zmax)
c
c         denver plotting system
c
20      continue
c20        if(iplotr.ge.9) go to 30
          if(ic.eq.2) go to 78
          ksw=2
c          if(iplotr.eq.1.or.iplotr.eq.4) go to 75
          nxl=(xdata(npts)-xdata(1))/xscale*2.+.5
          if(nxl.lt.2) nxl=2
          nyl=nxl
77        call scale3(xdata(1),xdata(npts),nxl,xminp,xmaxp,xdel)
          call scale3(zmin,zmax,nyl,yminp,ymaxp,ydel)
          xd(1)=xminp
          xd(2)=xmaxp
          yd(1)=yminp
          yd(2)=ymaxp
          xs(1)=(xmaxp-xminp)/xscale
          ys(1)=xs(1)*.5
          if(ys(1).gt.18.) ys(1)=18.
          go to 80
78        nyl=(zmax-zmin)/xscale*vexag*4.+.5
          if(nyl.lt.2) nyl=2                 
135       call scale3(zmin,zmax,nyl,yminp,ymaxp,ydel)
          yd(1)=ymaxp
          yd(2)=yminp
          xs(1)=(xmaxp-xminp)/xscale
          ys(1)=(ymaxp-yminp)/xscale*vexag   
80        call pltset(iplotr,xs(4),ys(4),jdim)
          yyss=ys(4)
          xs(2)=0.
          ys(2)=0.
          xs(3)=2.
          ys(3)=2.
          if((xs(1)+xs(3)+.2).gt.xs(4)) go to 85
          if((ys(1)+ys(3)+.2).gt.ys(4)) go to 85
          go to 87
85        xs(4)=amin1(xs(4),13.)
          if(ic.eq.1) go to 86
          ys(4)=(ymaxp-yminp)/(xmaxp-xminp)*xs(4)*vexag+1.2
          if(ys(4).gt.yyss) ys(4)=yyss
          go to 88
86        ys(4)=amin1(ys(4),10.)
88        xs(3)=1.0
          ys(3)=1.0
          xs(1)=xs(4)-xs(3)-.2
          ys(1)=ys(4)-ys(3)-.2
          go to 89
87        xs(4)=xs(1)+xs(3)+.2
          ys(4)=ys(1)+ys(3)+.2
89        if(revx.ne.'y') go to 90
          xd(1)=xmaxp
          xd(2)=xminp
90        call scale(xd,yd,xs,ys,4,icode)
          if(icode.ne.0) go to 99
          call neatl
          call xaxis(xd,yd,xs,xdel,2,.08,'(f10.2)',10)
          call yaxis(yd,xd,ys,ydel,2,.08,'(f10.2)',10)
        call vchar(2.,.3,xsym,40,3,.15,.0,.0,.0)
          if(ic.eq.2) return
c
c         input profile
c
        call vchar(2.,.1,sym,40,3,.15,0.,0.,0.)
          call line(xdata,zbuf,npts,0,0)
        jdim(2)=idim
69      call endpt(jdim)
          return
70      jdim(2)=idim
        call endpt(jdim)
        return
99        write(iprint,1200)
1200      format(' unable to scale plotter')
          stop
          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 minmax(a,n,amin,amax)
      dimension a(1)
      amin=a(1)
      amax=amin
      do 1 i=2,n
      amin=amin1(amin,a(i))
      amax=amax1(amax,a(i))
    1 continue
      return
      end
      subroutine scale3(xmin,xmax,n,xminp,xmaxp,dist)
c--see comm. acm, vol.16, no.10, oct 1973, pp.639-640.
      dimension vint(5)
      data vint(1),vint(2),vint(3),vint(4),vint(5)/1.,2.,5.,10.,20./
      if(xmin.lt.xmax.and.n.gt.1) go to 10
      write(6,99999)
99999 format(34h improper input supplied to scale3)
      return
   10 del=.00002
      fn=n
      a=(xmax-xmin)/fn
      al=alog10(a)
      nal=al
      if(a.lt.1.) nal=nal-1
      b=a/10.**nal
      do 20 i=1,3
      if(b.lt.(vint(i)+del)) go to 30
   20 continue
      i=4
   30 dist=vint(i)*10.**nal
      fm1=xmin/dist
      m1=fm1
      if(fm1.lt.0.) m1=m1-1
      if(abs(float(m1)+1.-fm1).lt.del) m1=m1+1
      xminp=dist*float(m1)
      fm2=xmax/dist
      m2=fm2+1.
      if(fm2.lt.(-1.)) m2=m2-1
      if(abs(fm2+1.-float(m2)).lt.del) m2=m2-1
      xmaxp=dist*float(m2)
      np=m2-m1
      if(np.le.n) go to 40
      i=i+1
      go to 30
c
c         next 3 statements comented out by r.godson,usgs,denver,co.
c
c   40 nx=(n-np)/2
c      xminp=xminp-float(nx)*dist
c      xmaxp=xminp+float(n)*dist
40      if(xminp.gt.xmin) xminp=xmin
      if(xmaxp.lt.xmax) xmaxp=xmax
      return
      end
          subroutine splin1(m,h,x,y,a,b,c,t,d,p,s)
c--one dimensional cubic spline interpolation
c
c        by  w.l.anderson, u.s. geological survey,
c             denver, colorado
c
c  parms--- m= number of data points .gt. 2
c           h= equal interval option when h.gt.0.
c              (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
c                 a,b,&c.
c                 also used as work arrays during execution.
c           t= type of boundary condition supplied in d array.
c              use t=1 if 1st derivatives given at end points, or
c              t=0 if 2nd derivatives given at end points.
c           d= boundary array (dim 2) at point 1 and m
c              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
          real   x(1),y(2),a(2),b(100),c(100),d(2),
     & p(100),s(100),mul
          integer  t
          if(t.lt.0.or.t.gt.1.or.m.lt.3) go to 999
          n=m-1
          if(t.eq.0) go to 20
c--1st derivative boundaries given
          ne=n-1
          if(h) 1,11,1
c--equal spacing h .gt. 0. and t=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         fa=1./b(1)
          c(1)=c(1)*fa
          p(1)=p(1)*fa
          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+t)=p(ne)
          i=ne-1
5         a(i+t)=p(i)-c(i)*a(i+t+1)
          i=i-1
           if(i.ge.1) go to 5
          if(t.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 t=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) 21,31,21
c--equal spacing h .gt. 0 and t=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 t=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
      subroutine namemc(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
        subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c       assigns values to proper variables
c       variables are passed to program acdepth through common blocks
c       numr=position in the array var where real varialbes start
c       numa=position in the array var where arrays start
c       nnvar=number of variables in program acdepth
c
        parameter (nnvar=20,numr=10)
        character*6 pvar,var(nnvar)
        character*56 tvar, kvar, cfmt
        character*1 crt(27000),revx,sym
        logical chv  
        dimension iprof(50),zbuf(3000),xdata(3000)
        common/mintc1/idplot, ind, revx, sym, crt
        common/mintco/zbuf, xdata, npts, iplotr, deltx, xscale, vexag,
     1  delx,dely
        common/parms/ichn, n, eps, iskip, incr, nskip, zmin, zmax, ipr,
     1  filnam, ifmt, name, iprof, outfile
        character*40 idplot,ind
        character*80 ifmt
        character*50 filnam, outfile
        character*4 name
        data var/'ichn', 'n', 'iskip', 'incr', 'nskip', 'ipr',
     1  'iplotr','iprof', 'npts', 'eps', 'zmin', 'zmax', 'xscale',
     2  'vexag', 'deltx', 'idplot', 'filnam', 'ifmt', 'name', 'outfil'/
        numa= 8
        do 190 i=1,nnvar
        if (pvar .ne. var(i)) goto 190
        if (.not. chv) then
c
c       non character value
c       right justify the number in variable kvar
c
        m = 57 - nn
        im = m - 1
        kvar(m:56) = tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write(cfmt, 50) im, nn
50      format ('(',i2,'x,i',i2,')')
        read(kvar, cfmt) jvar
        else
c
c       real value
c
        write (cfmt, 60) im, nn
60      format ('(',i2,'x,g',i2, '.0)')
        read (kvar, cfmt) xvar
        endif
        endif
        goto (101,102,103,104,105,106,107,108,109,110,111,112,113,
     1  114,115,116,117,118,119,120),i
101     ichn = jvar
        goto 200
102     n = jvar
        goto 200
103     iskip = jvar
        goto 200
104     incr = jvar
        goto 200
105     nskip = jvar
        goto 200
106     ipr = jvar
        goto 200
107     iplotr = jvar
        goto 200
108     iprof(inum)=jvar
        goto 200
109     npts = jvar
        goto 200
110     eps = xvar
        goto 200
111     zmin = xvar
        goto 200
112     zmax = xvar
        goto 200
113     xscale = xvar
        goto 200
114     vexag = xvar
        goto 200
115     deltx = xvar
        goto 200
116     idplot = tvar(1:nn)
        goto 200
117     filnam = tvar(1:nn)
        goto 200
118     ifmt = tvar(1:200)
        goto 200
119     name = tvar(1:nn)
        goto 200
120     outfil = tvar(1:nn)
        goto 200
190     continue
        write(*,*)' error in namelist - ',pvar,' variables not included'
        stop
200     nvar = i
        return
        end
