	subroutine pltdvr(ier)
c  plot driver for observed and calculated data
        common /field/  nf(10),f(5120)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1		iave,jave,vex,istatn,iparsh,labelv,numbod
	common /plot1/  iplotr,sizea,sizet,sizel,
     1		title1,title2,title3,
     1		xscale,xxx(2),adelx1,lintx,itx(5),
     1		zscale,zzz(2),adelz1,lintz,itz(5),
     1		gscale,ggg(2),adelg1,lintg,itg(5),
     1		hscale,hhh(2),adelh1,linth,ith(5)
	common /window/ xwind(2,4),ywind(2,4)
	common /vuport/ xview(4,4),yview(4,4)
	common /scal1/  xp(4),yp(4),dxp(2),dyp(2)
	common /pltdev/ jplotr,xboard,yboard,dum(2)
c
	dimension nlab(3),it1(14)
	character title1*56,title2*56,title3*56,titl*56
	logical grvswt,magswt
	equivalence (titl,it1)
        data ix,iobs,icalc,ioff/1,3,4,5/
c
c  initialize plot
        call clrscr
	if(grvswt) mgswt=0
	if(magswt) mgswt=1
	if(grvswt .and. magswt) mgswt=-1
	startx=0.01
	starty=0.01
	nplot=2
	if(mgswt.eq.-1) nplot=3
	call setwin(mgswt,xxx,zzz,ggg,hhh)
	call openvu(0,ier)
	call setvu(startx,starty,nlab,nlabel,ier)
	if(ier.gt.0) then
	  print *,'pltdvr:setvu error=',ier
c	  call waiter(2.0)
	  go to 99
	endif
c
c  plot model
	call openvu(1,ier)
	if(ier.eq.0) call pltmod
c
c  gravity plot
	if(grvswt) then
	iview=2
	np=nobs
	call openvu(iview,ier)
	if(ier.eq.0) then
	ltyp=0
	call pltgrf(f(nf(ix)),f(nf(icalc)),np,ltyp,iview)
c    character 6 is a '+'
	ltyp=-6
	call pltgrf(f(nf(ix)),f(nf(iobs)),np,ltyp,iview)
	endif
	endif
c
c  magnetic plot
	if(magswt) then
	iview=3
	np=mobs
	call openvu(iview,ier)
	if(ier.eq.0) then
	ltyp=0
	call pltgrf(f(nf(ix+ioff)),f(nf(icalc+ioff)),np,ltyp,iview)
	ltyp=-11
c     character 11 is a triangle
	call pltgrf(f(nf(ix+ioff)),f(nf(iobs+ioff)),np,ltyp,iview)
	endif
	endif
c
c  main titles
	if(ier.eq.0) then
	xc=xp(3)+0.5*xp(1)
	yt=yp(3)+yp(1)+((nlabel-1)+0.5)*(1.5*sizel)
	size=sizel
	do 10 i=1,3
	if(i.gt.1) size=.75*sizel
	if(nlab(i).gt.0) then
	  if(i.eq.1) titl(1:56)=title1(1:56)
	  if(i.eq.2) titl(1:56)=title2(1:56)
	  if(i.eq.3) titl(1:56)=title3(1:56)
	  xt=xc-0.5*float(nlab(i))*size
	  call vchar(xt,yt,it1,nlab(i),3,size,0.0,0.0,0.0)
	  yt=yt-1.5*sizel
	endif
10      continue 
	endif
c
c  close plot
99	call openvu(-1,ier)
        call clrscr
	return
	end
c******************************************************************************
	subroutine isodvr
c  plot the calculated curves for isolated groups of bodies	
        common /field/  nf(10),f(5120)
	common /swork/  a(8420)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1               iave,jave,vex,istatn,iparsh,labelv,numbod
	common /window/ xwind(2,4),ywind(2,4)
	dimension includ(20)
	logical grvswt,magswt
	data itty/5/
c
	m=nobs+mobs
c  put curves in the work array
10	print *,' enter number of curves to be plotted (max=7)'
	read(itty,*,end=999,err=999) ncurv
	if(ncurv.lt.1 .or. ncurv.gt.7) go to 999
	do 20 i=1,m*ncurv
	  a(i)=0.0
20      continue
	ig=1
	ih=nobs+1
	do 40 icurv=1,ncurv
	i1=icurv
	print *,' enter nbody and body indices for curve', icurv
	read(itty,*,end=10,err=10) n,(includ(i),i=1,n)
	  do 30 i=1,n
	  if(grvswt) then
	    call sumfwd(0,includ(i),nobs,a(ig))
	    call ammi(nobs,a(ig),ywind(1,2),ywind(2,2),mn,mx,i1)
	  endif
	  if(magswt) then
	    call sumfwd(1,includ(i),mobs,a(ih))
	    call ammi(mobs,a(ih),ywind(1,3),ywind(2,3),mn,mx,i1)
	  endif
30        continue 
	ig=ig+m
	ih=ih+m
40      continue 
c
	call plota(1,nobs,1,mobs,m,ncurv,a)
999	return
	end
c******************************************************************************
	subroutine sumfwd(mgswt,ibody,np,pf)
        common /field/  nf(10),f(5120)
	common /parm1/  nvert,nbody,nobs,mobs,dum(9)
	common /model/  xb(200),zb(200),rho(50),sus(50),pl(50,2),
     1               rem(50,3),dum1(150)
	common /vertx/  icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	dimension xb1(200),zb1(200),pf(np)
	do 10 i=1,ncrn(ibody)
	  xb1(i)=xb(icrn(loca(ibody)+i))
	  zb1(i)=zb(icrn(loca(ibody)+i))
10      continue 
	if(mgswt.eq.0) then
	  rs=rho(ibody)
	  ix=1
	  iz=ix+nobs
	  else
	  rs=sus(ibody)	
	  ix=5*nobs+1
	  iz=ix+mobs
	endif
	call pickmg(mgswt,ibody,rs,ncrn(ibody),xb1,zb1,
     1           np,f(ix),f(iz),pf)
	return
	end
c******************************************************************************
	subroutine advr(igs,ngrv,ihs,nmag,no,npart,a,e)
c  driver for the partial derivative and error vector plots
        common /field/  nf(10),f(5120)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1               iave,jave,vex,idepth,iparsh,indiv,numbod
	common /freel/  ifill(4),nx,nz,nr,ns,ix(20),iz(20),ir(20),is(20)
	common /window/ xwind(2,4),ywind(2,4)
	dimension xx(2),zz(2),gg(2),hh(2)
	dimension a(no,npart),e(no)
	logical grvswt,magswt
c
	mgswt=0
	if(magswt) mgswt=1
	if(grvswt .and. magswt) mgswt=-1
c  set x window
	xx(1)= 1.0e5
	xx(2)=-1.0e5
	if(grvswt) then
	  igx=nf(1)+(igs-1)
	  xx(1)=f(igx)
	  xx(2)=f(igx+ngrv-1)
	endif
	if(magswt) then
	  ihx=nf(6)+(ihs-1)
	  xx(1)=amin1 (f(ihx),xx(1))
	  xx(2)=amax1 (f(ihx+nmag-1),xx(2))	
	endif
c
c  x partials
	if(nx.gt.0) then
	do 10 j=1,nx
	j1=j
	if(grvswt) call ammi(ngrv,a(1,j),gg(1),gg(2),
     1          mn,mx,j1)
	if(magswt) call ammi(nmag,a(ngrv+1,j),hh(1),hh(2),
     1          mn,mx,j1) 
10      continue 
	call setwin(mgswt,xx,zz,gg,hh)
	call plota(igs,ngrv,ihs,nmag,no,nx,a)
	endif
c
c  z partials
	if(nz.gt.0) then
	js=nx+1
	je=nx+nz
	do 20 j=js,je
	j1=j
	if(grvswt) call ammi(ngrv,a(1,j),gg(1),gg(2),
     1          mn,mx,j1)
	if(magswt) call ammi(nmag,a(ngrv+1,j),hh(1),hh(2),
     1          mn,mx,j1) 
20      continue 
	call setwin(mgswt,xx,zz,gg,hh)
	call plota(igs,ngrv,ihs,nmag,no,nz,a(1,js))
	endif
c
c  density and susceptibility
	if(nr+ns.gt.0) then
	js=nx+nz+1
	do 30 j=js,nx+nz+nr
	j1=j
	call ammi(ngrv,a(1,j),     gg(1),gg(2),mn,mx,j1)
30      continue 
	do 40 j=js+nr,npart
	call ammi(nmag,a(ngrv+1,j),hh(1),hh(2),mn,mx,j1) 
40      continue
	call setwin(mgswt,xx,zz,gg,hh)
	nrs=nr+ns
	call plota(igs,ngrv,ihs,nmag,no,nrs,a(1,js))
	endif
c
c  error vector
	if(grvswt) call ammi(ngrv,e(1),     gg(1),gg(2),mn,mx,1)
	if(magswt) call ammi(nmag,e(ngrv+1),hh(1),hh(2),mn,mx,1)
	call setwin(mgswt,xx,zz,gg,hh)
	call plota(igs,ngrv,ihs,nmag,no,1,e)
	return
	end
c******************************************************************************
	subroutine plota(igs,ngrv,ihs,nmag,m,ncurv,a)
c  plot multiple curves located in 'a'.
c  the functional values in 'a' are column-wise arranged grv then mag
c  with location i=1 corresponding to igx or ihx.
c  igs,ihs refer to the starting locations in field common.
        common /field/  nf(10),f(5120)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1               iave,jave,vex,idepth,iparsh,indiv,numbod
	common ifill(4),nx,nz,nr,ns,ix(20),iz(20),ir(20),is(20)
	dimension a(m,ncurv)
	logical grvswt,magswt
c
	if(ngrv+nmag.ne.m) return
	igx=nf(1)+(igs-1)
	ihx=nf(6)+(ihs-1)
        call clrscr
	call openvu(0,ier)
	if(ier.eq.0) call setvu(startx,starty,nlab,nlabel,ier)
	if(ier.eq.0) call openvu(1,ier)
	if(ier.gt.1) go to 88
	call pltmod
	if(grvswt) then
	  call openvu(2,ier)
	  if(ier.gt.0) go to 88
	  do 10 j=1,ncurv
	  ltyp=j-1
	  call pltgrf(f(igx),a(1,j),ngrv,ltyp,2)
10        continue 
	endif
	if(magswt) then
	  call openvu(3,ier)
	  if(ier.gt.0) go to 88
	  ih=ngrv+1
	  do 20 j=1,ncurv
	  ltyp=j-1
	  call pltgrf(f(ihx),a(ih,j),nmag,ltyp,3)
20        continue 
	endif
88	call openvu(-1,ier)
        call clrscr
99	return
	end
c******************************************************************************
	subroutine pltmod
        common /field/  nf(10),f(5120)
	common /model/  xb(200),zb(200),fill(500)
	common /vertx/  icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1		iave,jave,vex,istatn,iparsh,labelv,numbod
	common /plot1/  iplotr,sizea,sizet,sizel,
     1		it1(14),it2(14),it3(14),
     1		xscale,xxx(2),adelx,lintx,itx(5),
     1		zscale,zzz(2),adelz,lintz,itz(5),
     1		gscale,ggg(2),adelg,lintg,itg(5),
     1		hscale,hhh(2),adelh,linth,ith(5)
	common /scal1/  xp(4),yp(4),dxp(2),dyp(2)
c
	dimension zlmt(2),ify(5),xb1(100),zb1(100)
	logical grvswt,magswt
	character fmty*20,titlz*20,body*3
	equivalence (fmty,ify),(titlz,itz)
c
c  draw model cross-section
	xsca=xp(1)/(dxp(2)-dxp(1))
	ysca=yp(1)/(dyp(2)-dyp(1))
	szlabl=1.25*sizea
	sizekm=sizea/xsca
	do 30 j=1,nbody
	do 10  i=1,ncrn(j)
	i2=icrn(loca(j)+i)
	xb1(i)=xb(i2)
10	zb1(i)=zb(i2)
	i1=ncrn(j)+1
	xb1(i1)=xb1(1)
	zb1(i1)=zb1(1)
	call line(xb1,zb1,i1,0,200)
c
c  number bodies
	if(numbod.eq.1) then
	  write(body,20) j
20	  format(i2)
	  zlmt(1)=dyp(2)
	  zlmt(2)=dyp(1)
	  nch=leftj(body)
	  call inpoly(a,b,ncrn(j),xb1,zb1,dxp,zlmt,sizekm,nch,ier)
	  if(ier.eq.0 .and. sizea.gt.0.0) then
c           convert location to plotter inches
	    a=xp(3)+xsca*(a-dxp(1))
	    b=yp(3)+ysca*(b-dyp(1))
	    call text(a,b,body,szlabl,0.0)
	    endif
	endif
30	continue
c
c  plot observation locations
	if(istatn.eq.1) then
	  if(grvswt) then
	  ich=6
	  call symbol(f(nf(1)),f(nf(2)),nobs,ich,sizea,0.0)
	  endif
	  if(magswt) then
	  ich=11
	  call symbol(f(nf(6)),f(nf(7)),mobs,ich,sizea,0.0)
	  endif
	endif
c
c  depth axis annotation
	adely=adelz
	call setax(dyp,adely,15,nchary,fmty)
	call yaxis(dyp,dxp,yp,adely,lintz,sizea,ify,nchary)	
c        if(iplotr.ne.1) then
	  nch=leftj(titlz)
	  yt=0.5*yp(1)+yp(3)-float(nch/2)*sizet
	  xt=xp(3)-nchary*sizea-sizet
	  call text(xt,yt,titlz,sizet,90.0)
c        endif
c
c  label vertices with a diamond and number 
	if(labelv.eq.1) then
	ich=2
	tqsz=0.75*sizea
	hsz =0.5 *sizea
	  do 40 iv=1,nvert
	  call symbol(xb(iv),zb(iv),1,ich,hsz,0.0)
	  a=xp(3)+xsca*(xb(iv)-dxp(1))
	  b=yp(3)+ysca*(zb(iv)-dyp(1))
         write(body,35) iv
35       format(i3)
	  b=b-tqsz
	  call text(a,b,body,tqsz,0.0)
40    continue
	endif
	return
	end
c******************************************************************************
	subroutine inpoly(a,b,n,x,y,xlmt,ylmt,size,nsym,ier)
c  place labeling box inside a polygon, clip at xlmt, ylmt
	dimension xlmt(2),ylmt(2),x(n),y(n)
	logical inside
	ier=1
	ntst=20
	hsize = size/2.0
	xsize = (nsym+1)*(5.0/7.0)*size
	call ammi(n,x,xmn,xmx,mn,mx,0)
	call ammi(n,y,ymn,ymx,mn,mx,0)
	if(xlmt(1).gt.xmn) xmn=xlmt(1)
	if(xlmt(2).lt.xmx) xmx=xlmt(2)
	if(ylmt(1).gt.ymn) ymn=ylmt(1)
	if(ylmt(2).lt.ymx) ymx=ylmt(2)
	if(xmx-xmn.le.0.0) return
	if(ymx-ymn.le.0.0) return
	as = (xmn+xmx)/2.0
	bs = (ymn+ymx)/2.0
	a=as
	b=bs
c
	dx=(xmx-xmn)/float(ntst-1)
	dist=0.0
	b1=bs-hsize
	b2=b1+size
	do 10 itst=0,ntst
	dist=-sign(1.0,dist)*(float(itst)*dx)
	a=as+dist
	a1=a-hsize
	a2=a1+xsize
	if(inside(n,x,y,a1,b1) .and. inside(n,x,y,a2,b2)) go to 99
10    continue
c
	dy=(ymx-ymn)/float(ntst-1)
	dist=0.0
	a1=as-hsize
	a2=a1+xsize
	do 20 itst=0,ntst
	dist=-sign(1.0,dist)*(float(itst)*dy)
	b=bs+dist	
	b1=b-hsize
	b2=b1+size
	if(inside(n,x,y,a1,b1) .and. inside(n,x,y,a2,b2)) go to 99
20    continue
	return
99	ier=0
	return
	end
c******************************************************************************
	subroutine pltgrf(x,y,np,ltyp,iview)
c  line print = 0 annotates and draws the first curve,
c            > 0 draws curves with various dashing patterns
c            < 0 draws symbols
	common /parm1/  nvert,nbody,ifill(2),grvswt,magswt,
     1		iave,jave,vex,istatn,iparsh,labelv,numbod
	common /plot1/  iplotr,sizea,sizet,sizel,
     1		it1(14),it2(14),it3(14),
     1		xscale,xxx(2),adelx,lintx,titlx,
     1		zscale,zzz(2),adelz,lintz,titlz,
     1		gscale,ggg(2),adelg,lintg,titlg,
     1		hscale,hhh(2),adelh,linth,titlh
	common /scal1/  xp(4),yp(4),dxp(2),dyp(2)
	dimension ifx(5),ify(5)
	dimension x(np),y(np)
	character*20 fmtx,fmty,titlx,titlz,titlg,titlh,tx,ty
	logical grvswt,magswt
	equivalence (fmtx,ifx),(fmty,ify)
c
	if(ltyp) 300,100,200
100	adelx1=adelx
	adely=adelg
	linty=lintg
	tx(1:20)=titlx(1:20)
	if(iview.eq.3 .and. (grvswt.and.magswt)) tx=' '
	ty(1:20)=titlg(1:20)
	if(iview.eq.3) then
	adely=adelh
	linty=linth
	ty(1:20)=titlh(1:20)
	endif
c
	call setax(dyp,adely, 30,nchary,fmty)
	call yaxis(dyp,dxp,yp,adely, linty,sizea, ify,nchary)
	nchy=leftj(ty)
	if(nchy.gt.0) then
	  xt=xp(3)-nchary*sizea-sizet
	  yt=yp(3)+0.5*yp(1)-(nchy/2)*sizet
	  call text(xt,yt,ty,sizet,90.0)
	endif
c
	call setax(dxp,adelx1,50,ncharx,fmtx)
	if(iview.eq.3 .and. grvswt) ncharx=0
	call xaxis(dxp,dyp,xp,adelx1,lintx,sizea,ifx,ncharx)
	nchx=leftj(tx)
	if(nchx.gt.0) then
	  xt2=xp(3) + 0.5*xp(1) - (nchx/2)*sizet
	  yt2=yp(3)-ncharx*sizea-sizet
	  call text(xt2,yt2,tx,sizet,0.0)
	endif
c
c  plot calculated function
200	ltyp1=mod(ltyp,8)
	call line(x,y,np,0,ltyp1+300)
	return
c
c  plot observed data
300	ich=iabs(ltyp)
	call symbol(x,y,np,ich,sizea,0.0)
	return
	end
c******************************************************************************
	subroutine setwin(mgswt,xxx,zzz,ggg,hhh)
c  assign window limits either from data or variables xxx,zzz,ggg,hhh.
c  the triple letter input variables have priority.
c  mgswt=-1 both fields, =0 gravity only, =1 magnetics only.
        common /field/  nf(10),f(5120)
	common /model/  xbody(200),ybody(200),fill(500)
	common /vertx/  fill2(800),nbord,ibord(200)
	common /parm1/  ifill1(2),nobs,mobs,ifill2(5),istatn,ifill3(3)
	common /window/ xwind(2,4),ywind(2,4)
	dimension xxx(2),zzz(2),ggg(2),hhh(2)
	data ix,iz,iobs,icalc,ioff /1,2,3,4,5/
c
	do 10 j=1,3
	  xwind(1,j)= 1.e10
	  xwind(2,j)=-1.e10
	  ywind(1,j)= 1.e10
	  ywind(2,j)=-1.e10
10    continue
c
c  x windows
	if(mgswt.lt.1 .and. nobs.gt.0)
     1  call ammi(nobs,f(nf(ix)),     xwind(1,2),xwind(2,2),mn,mx,1)
	if(mgswt.ne.0 .and. mobs.gt.0)
     1  call ammi(mobs,f(nf(ix+ioff)),xwind(1,3),xwind(2,3),mn,mx,1)
c
	if((xxx(1).ne.0.0 .or. xxx(2).ne.0.0) .and.
     1	xxx(1).lt.xxx(2)) then
	  xwind(1,1)=xxx(1)
	  xwind(2,1)=xxx(2)
	else
	  xwind(1,1)=amin1(xwind(1,2), xwind(1,3))
	  xwind(2,1)=amax1(xwind(2,2), xwind(2,3))
	endif
c
	do 20 j=2,3
	  xwind(1,j)=xwind(1,1)
	  xwind(2,j)=xwind(2,1)
20    continue
c
c  y windows
c     gravity
	if(mgswt.lt.1 .and. nobs.gt.0) then
	  if(istatn.eq.1)  call ammi(nobs,f(nf(iz)),
     1                 ywind(1,1),ywind(2,1),mn,mx,2)
	  call ammi(nobs,f(nf(iobs)),  ywind(1,2),ywind(2,2),mn,mx,2)
	  call ammi(nobs,f(nf(icalc)), ywind(1,2),ywind(2,2),mn,mx,2)
	endif
c    magnetics
	if(mgswt.ne.0 .and. mobs.gt.0) then
	  if(istatn.eq.1) call ammi(mobs,f(nf(iz+ioff)),
     1                 ywind(1,1),ywind(2,1),mn,mx,2)
	  call ammi(mobs,f(nf(iobs+ioff)),  ywind(1,3),ywind(2,3),mn,mx,2)
	  call ammi(mobs,f(nf(icalc+ioff)), ywind(1,3),ywind(2,3),mn,mx,2)
	endif
c   model
	do 30 k=1,nbord	
	  yb=ybody(ibord(k))
	  if(yb.lt.ywind(1,1)) ywind(1,1)=yb	
	  if(yb.gt.ywind(2,1)) ywind(2,1)=yb
30    continue
c   minimum range
	do 40 j=2,3
	  yrange=abs(ywind(2,j)-ywind(1,j))
	  if(yrange.eq.0.0) yrange=1.0
	  ywind(2,j)=ywind(2,j)+.05*yrange
	  ywind(1,j)=ywind(1,j)-.05*yrange
40    continue
c
c  specified limits
	if((zzz(1).ne.0.0 .or. zzz(2).ne.0.0)) then
	  ywind(1,1)=zzz(1)
	  ywind(2,1)=zzz(2)
	endif
	if((ggg(1).ne.0.0 .or. ggg(2).ne.0.0) .and.
     1	ggg(1).lt.ggg(2)) then
	  ywind(1,2)=ggg(1)
	  ywind(2,2)=ggg(2)
	endif
	if((hhh(1).ne.0.0 .or. hhh(2).ne.0.0) .and.
     1	hhh(1).lt.hhh(2)) then
	  ywind(1,3)=hhh(1)
	  ywind(2,3)=hhh(2)
	endif
c
c  ad hoc... z positive down in plot area
	if(ywind(1,1).lt.ywind(2,1)) then
	  tmp=ywind(1,1)
	  ywind(1,1)=ywind(2,1)
	  ywind(2,1)=tmp
	endif
	return
	end
c******************************************************************************
	subroutine setvu(startx,starty,nlab,nlabel,ier)
c  Set dimensions for the cross-section, gravity, and magnetic 
c  viewports.  When the scaling defaults, the associated viewport is
c  a set fraction of the total device width.  Viewport parameters in 
c  plotter units: 1,2 are minimum,maximum coordinates
c  3,4 are left,right interior margin widths for a scaled data area.
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,
     1		iave,jave,vex,istatn,iparsh,labelv,numbod
	common /plot1/  iplotr,sizea,sizet,sizel,
     1		title1,title2,title3,
     1		xscale,xxx(2),adelx,lintx,itx(5),
     1		zscale,zzz(2),adelz,lintz,itz(5),
     1		gscale,ggg(2),adelg,lintg,itg(5),
     1		hscale,hhh(2),adelh,linth,ith(5)
	common /pltdev/ jplotr,xboard,yboard,dum(2)
	common /window/ xwind(2,4),ywind(2,4)
	common /vuport/ xview(4,4),yview(4,4)
	character title1*56,title2*56,title3*56
	logical grvswt,magswt
	dimension scale(3),nlab(3)
c
	nport=3
	nplot=2
	if(grvswt .and. magswt) nplot=3
	ier=1
	if(xboard.le.0.0 .or. yboard.le.0.0) return
	do 10 i=1,nport
	  scale(i)=0.0
10    continue
	xsca=0.0
	if(xscale.gt.0.0) xsca=    xscale
	if(zscale.gt.0.0) scale(1)=zscale
	if(gscale.gt.0.0) scale(2)=gscale
	if(hscale.gt.0.0) scale(3)=hscale
c
c  labelling space
c        sizet1=sizet
c        sizel1=sizel
c        nchar=13
        nlabel=0
c        if(iplotr.eq.1) then
c   conserve space on video terminals
	  sizet1=0.0
	  sizel1=0.0
	  nchar=7
c        else
	  nlab(1)=leftj(title1)
	  nlab(2)=leftj(title2)
	  nlab(3)=leftj(title3)
	  do 20 i=1,3
	  if(nlab(i).gt.0) nlabel=nlabel+1
20    continue
c        endif
        xymarg=3.0*sizet1+nchar*sizea
	if(xymarg.lt.0.0) xymarg=0.0
	startx=0.02*xboard
	starty=0.02*yboard
c
c  x dimensions
	if(2.0*startx.gt.xboard) startx=0.0
	do 30 j=1,3
	xview(1,j)=startx
	xview(2,j)=xboard-startx
	xview(3,j)=xymarg
	xview(4,j)=sizea
	do 30 i=1,4
30	yview(i,j)=0.0
c
c  rightside margins with scaling included
	wdata=abs(xwind(2,1)-xwind(1,1))
	if(xsca.gt.0.0) then
	  xp1 = wdata/xsca
	  xview(4,1) = xview(2,1)-(xview(1,1)+xview(3,1)+xp1)
	  if(xview(4,1).lt.0.0) then
	    xview(4,1)=sizea
	    xwind(2,1) = xsca*p1(xview(1,1)) + xwind(1,1)
	  endif
	else
	  ier=2
	  if(p1(xview(1,1)).le.0.0) return
	  if(xsca.eq.0.0) xsca=wdata/p1(xview(1,1))
	  ier=3
	  if(xsca.le.0.0) return
	endif
c
	do 40 j=2,3
	xwind(2,j)=xwind(2,1)
	do 40 i=1,4
40	xview(i,j)=xview(i,1)
	if(scale(1).le.0.0) scale(1)=xsca
	if(vex.le.0.0) vex=1.0
	scale(1)=scale(1)/vex
c
c  model cross-section dimensions
c
	ydef=(yboard-2.*starty)/float(nport)
        ysep=0.1*ydef
	wdata=abs(ywind(2,1)-ywind(1,1))
c
c  viewport size
	yview(1,1)=starty
	yview(2,1)=ydef
	yview(3,1)=0.0
	yview(4,1)=ysep
	if(xscale.gt.0.0 .or. zscale.gt.0.0) then
	  yview(2,1) = yview(1,1)+ysep+wdata/scale(1)
	  ylimit = 0.7*(yboard-2.0*starty)
	  if(yview(2,1)-yview(1,1).gt.ylimit) yview(2,1)=ylimit
	endif
c
c  cutoff bottom of the window to maintain scaling
	yview(3,1)=yview(2,1)-yview(1,1)-yview(4,1)-wdata/scale(1)
	if(yview(3,1).lt.0.0) then
	  yview(3,1)=0.0
	  ywind(1,1)=ywind(2,1)+scale(1)*p1(yview(1,1))
	endif
c
c  graph dimensions, port 2 is gravity, port 3 is magnetics
c
	ydef=((yboard-starty)-yview(2,1))/(nplot-1)
	do 50 j=2,nport
	yview(1,j)=yview(2,j-1)
	yview(2,j)=yview(1,j)
	if(j.eq.2 .and. .not.grvswt) go to 50
	if(j.eq.3 .and. .not.magswt) go to 50
c
	if(j.eq.3 .and. grvswt) xymarg=0.0
	if(j.eq.3 .or. (j.eq.2 .and. .not.magswt))
     1  ysep=nlabel*(1.5*sizel)
	yview(3,j)=xymarg
	yview(4,j)=ysep
c
c  viewport size
	wdata=abs(ywind(2,j)-ywind(1,j))
	if(scale(j).gt.0.0) then
	  yview(2,j)=yview(1,j)+yview(3,j)+yview(4,j)+wdata/scale(j)
	  if(yview(2,j).gt.yboard-starty) yview(2,j)=yboard-starty
	  else
	  yview(2,j)=yview(1,j)+ydef
	endif
c
	ier=4
	yp1=p1(yview(1,j))
	if(yp1.le.0.0) then
	  yview(3,j)=0.0
	  yview(4,j)=0.0	
	  yp1=p1(yview(1,j))
	endif
	if(yp1.lt.0.0) return
c  adjust window if necessary to maintain specified scaling
	if(scale(j).gt.0.0 .and. yp1.ne.wdata/scale(j))
     1  ywind(1,j)=ywind(2,j)-scale(j)*yp1
c
50	continue
	ier=0
	return
	end
c******************************************************************************
	function p1(a)
	dimension a(4)
	p1=a(2)-a(1)-a(3)-a(4)
	return
	end
c******************************************************************************
	subroutine openvu(iview,ier)
c  translate window and viewport info into the plot system form
c  and perform a scaling call to activate the viewport.
c
c  short summary of the plot system by G.I. Evenden and R.R. Wahl USGS.
c
c    >>pltset(jplotr,xboard,yboard,1) initializes device in inches.
c  x & y board are returned plotter dimensions.
c
c    >>scale(dxp,dyp,xp,yp,nopts,ierror)
c  Definition of plot area: xp(1:4),yp(1:4) units are inches:
c  xp(1)= data area width, 2=0 is linear scale, 3=left margin, 
c  4=plotter width.
c  Definition of data area: dxp(1:2),dyp(1:2) in data units:
c  dxp(1)=minimum x coordinate, dxp(2)=maximum x coordinate.
c  Vectors are clipped at the data area boundary.
c
c    >>line(x,y,n,icon,ipen) where x,y are data unit arrays n long, 
c  icon.ne.0 starts new line, and ipen selects dashing pattern.
c    >>xaxis,yaxis axis tic routines
c    >>endpt closes out the system
c
c  viewport coordinates in inches
c  window coordinates in data units (kilometers, mgal, etc)
c
	common /pltdev/ jplotr,xboard,yboard,xlimit,ylimit
	common /window/ xwind(2,4),ywind(2,4)
	common /vuport/ xview(4,4),yview(4,4)
	common /scal1/  xvu(4),yvu(4),xdata(2),ydata(2)
        dimension xp(4),yp(4),ie(2)
c
	if(iview.lt.0) then
          ie(2)=2
          call endpt(ie)
	  xboard=0.0
	  yboard=0.0
	  return
	endif
c
	if(iview.eq.0) then
c  initialize plot system, return device dimensions
	  call pltset(jplotr,xboard,yboard,1)
c  normally limit electrostatic plotter size to one paper width.
	  if(jplotr.eq.6 .and. 
     1       (xlimit.gt.xboard .or. ylimit.gt.yboard)) then
	    xlimit=72.
	    ylimit=20.
	  endif
	  if(xlimit.lt.xboard) xboard=xlimit
	  if(ylimit.lt.yboard) yboard=ylimit
	  return
	endif
c
c  open viewport
	call getvu(xview(1,iview),yview(1,iview),xp,yp,ier)
	if(ier.ne.0) print *,' openvu:getvu error number', ier
	if(ier.gt.0) return
	ier=0
	call scale(xwind(1,iview),ywind(1,iview),xp,yp,4,ier)
	if(ier.ne.0) then
	  print *,'openvu: scaling error'
	  return
	endif
	do 5 i=1,4
	xvu(i)=xp(i)
5	yvu(i)=yp(i)
	do 10 i=1,2
	xdata(i)=xwind(i,iview)
10	ydata(i)=ywind(i,iview)
	return
	end
c******************************************************************************
	subroutine getvu(xview,yview,xp,yp,ier)
c  input viewport data, returns 4 parameter data area specification
c  all values are positive and in plot units
	common /pltdev/ jplotr,xboard,yboard,dum(2)
	dimension xview(4),yview(4),xp(4),yp(4)
	tol=1.e-4
	xp(2)=0.0
	yp(2)=0.0
	ier=10000
	xp(4)=xboard
	yp(4)=yboard
	if(xp(4).le.0.0 .or. yp(4).le.0.0) return	
	ier=1
	xwidth = xview(2)-xview(1)
	if(xwidth.lt.0.0) go to 50
	ier=2
	xp(3) = xview(1) + xview(3)
	if(xp(3).lt.0.0) return
	ier=3
	xp(1) = xwidth - (xview(3)+xview(4))
	if(xp(1).lt.0.0) go to 50
	ier=0
c
50	ier=ier+10
	ywidth = yview(2)-yview(1)
	if(ywidth.lt.0.0) go to 99
	ier=ier+20
	yp(3) = yview(1) + yview(3)
	if(yp(3).lt.0.0) go to 99
	yp(1) = ywidth - (yview(3)+yview(4))
	ier=ier+30
	if(yp(1).lt.0.0) go to 99
	ier=0
c
99	if(xp(3).gt.tol) xp(3)=xp(3)-tol
	if(yp(3).gt.tol) yp(3)=yp(3)-tol
	if(xp(1)+xp(3).gt.xp(4) .and. xp(1).gt.tol) xp(1)=xp(1)-tol
	if(yp(1)+yp(3).gt.yp(4) .and. yp(1).gt.tol) yp(1)=yp(1)-tol
	xptest=xp(4)-xp(1)-xp(3)
	if(xptest.lt.0.0) then
	  xp(1)=xp(4)
	  xp(3)=0.0
	  ier=ier+100
	endif
	yptest=yp(4)-yp(1)-yp(3)
	if(yptest.lt.0.0) then
	  yp(1)=yp(4)
	  yp(3)=0.0
	  ier=ier+1000
	endif
	return
	end
c******************************************************************************
	subroutine setax(x,dx,maxint,nch,fmt)
c  adjust interval and labeling format
c  input x-min max range, optional dx-spacing and maxint-intervals
c  returns dx, nch-number of significant figures, fmt-labeling format
	dimension x(2)
	character fmt*(*)
	ixpn(r)=int(alog10(abs(r))+100.)-100
	fmt='(1pe13.6)'
	nch=13
	if(maxint.le.0) maxint=20
	dt=abs(x(2)-x(1))/float(maxint)
	if(dt.lt.1.e-20 .and. dt.gt.-1.e-20) go to 20
	if(dx.ne.0.) go to 5
	p10=sign(1.0,dt)*10.**ixpn(dt)
	t1=dt/p10
	if(t1.le.1.0) dx=p10
	if(t1.gt.1.0) dx=2.*p10
	if(t1.gt.2.0) dx=5.*p10
	if(t1.gt.5.0) dx=10.*p10
5	idecm=0
	m10=ixpn(dx)
	if(m10.lt.0) idecm=iabs(m10)
	iw=4+idecm
	m20=ixpn( amax1(abs(x(1)),abs(x(2))) )
	if(m20.ge.1) iw=iw+m20
	if(iw.gt.9 .or. idecm.gt.9) go to 11
	write(fmt,10) iw,idecm
10	format('(f',i1,'.',i1,')' )
	nch=iw
	if(idecm.eq.0) nch=nch-1
11	return
20	if(dx.eq.0.) dx=1.e20
	return
	end
c******************************************************************************
	subroutine adnois(nf,f,a,grvswt,magswt)
	dimension nf(10),a(400),f(1)
	character ifun*20
	logical grvswt,magswt
	nobs=nf(2)-nf(1)
	mobs=nf(7)-nf(6)
c
340	print  341
341	format(' add gaussian noise to calculated or observed :'$)
	read(5,37,end=999) ifun
37	format(a20)
	call cvc(ifun,2)
	iswt=-1
	if(ifun(1:2).eq.'ob') iswt=0
	if(ifun(1:2).eq.'ca') iswt=1
	if(iswt.eq.-1) go to 349
3400	print  342
342	format(' enter percent error :'$)
	read(5,*,end=3400,err=3400) perc
	print 343
343	format(' correlate the noise ?'$)
        icorr=noyes()
	if(icorr.eq.1) then
345	print 346
346	format(' enter correlation width (samples) :'$)
	read(5,*,end=345,err=345) iw
	if(iw.lt.2) iw=2
	endif
	noff=0
	n=nobs
	if(grvswt) go to 351
350	noff=5
	n=mobs
351	iadr=nf(3+iswt+noff)
	call ammi(n,f(iadr),fmn,fmx,mn,mx,1)
	sd=0.5*(perc/100.)*(fmx-fmn)
	if(sd.lt.1.e-20) then
	print *,' zero range function'
	go to 349
	endif
	do 354 i=1,n+iw
        call gaussd(sd,0.,a(i))
354	continue
	if(icorr.eq.1) then
	t=0.0
	do 360 i=1,iw
360	t=t+a(i)
	do 370 i=1,n
	ave=t/float(iw)
	t=t-a(i)+a(i+iw)	
370	a(i)=ave
	endif
	iadr1=iadr-1
	do 380 i=1,n
380	f(iadr1+i)=f(iadr1+i)+a(i)
	if(magswt .and. noff.eq.0) go to 350
	return
349	print  348
348	format(' no change')
999	return
	end
c******************************************************************************
        subroutine gaussd(sd,xmean,v)
c	save iseed
c	data is/1111111111/
c	if(iseed.eq.0) iseed=is
       data ix/0/,iy/0/,iz/0/
	a=0.
	do 1 i=1,12
1	a=a+ran3(ix,iy,iz)
	v=(a-6.0)*sd+xmean
	return
	end
c******************************************************************************
	subroutine trend(x,y,n,dc,dy)
	dimension x(n),y(n)
	dy1=(y(n)-y(1))/(x(n)-x(1))
	dc1=y(1)
	print  10,dc1,dy1
10	format(' dc and slope to zero ends:',2f10.4)
	call fit1(x,y,n,0.,dy2,dc2)
	print  11,dc2,dy2
11	format(' dc and slope to fit least square :',2f10.4)
	print  12
12	format(' enter desired dc and slope for removal ')
	read(5,*,end=999,err=999) dc,dy
	do 20 i=1,n
20	y(i)=y(i)-dc-(x(i)-x(1))*dy
	return
999	print  98
98	format(' no change')
	return
	end
c******************************************************************************
	subroutine drpsur(xo,dx,nobs,zl, xf,zf)
	common /vertx/ ifill(800), nbord, ibord(200)
	common /model/ xb(200),zb(200),dum(500)
	dimension xf(nobs),zf(nobs)
	if(nbord.le.0) go to 999
100	print *,' enter  start, ending vertices '
	print 1
1	format(' defining the top surface: '$)
	read(5,*,end=100,err=100) istart,iend
	zl = abs(zl)
c
c  js,je are ibord indices corresponding to vertices istart,iend. 
	js = 0
	je = 0
	do 10 i = 1, nbord
	  if (ibord(i).eq.istart) js=i
	  if (ibord(i).eq.iend)   je=i
10    continue
c
	do 20 i = js, je-1
	if(xb(ibord(i+1)).lt.xb(ibord(i))) then
	  print *,' top surface needs to be increasing in x'
	  return
	endif
20    continue
c
	do 40 i = 1, nobs
	xf(i) = xo + dx*float(i-1)
	zf(i) = 0.0
	if       (xf(i).lt.xb(ibord(js))) then
	  zf(i) = zb(ibord(js)) - zl
	  elseif (xf(i).gt.xb(ibord(je))) then
	  zf(i) = zb(ibord(je)) - zl
	endif
	if (xf(i).ge.xb(ibord(js)) .and. xf(i).le.xb(ibord(je))) then
	j = js
30	if (j.ge.je) go to 40
	k  = ibord(j)
	k1 = ibord(j+1)
	if (xf(i).ge.xb(k) .and. xf(i).lt.xb(k1)) then
	  xd = xb(k1) - xb(k)
	  if (xd.le.1.e-20) go to 40
	  zd = zb(k1) - zb(k)
	  x1 =  xf(i) - xb(k)
	  zf(i) = (x1*zd)/xd + (zb(k)-zl)
	else
	  j=j+1
	  go to 30
	endif
	endif
40	continue
c
	return
999	print  998
998	format(' no boundary array available, cannot drape')
	return
	end
c******************************************************************************
	subroutine svd(a,m,n,s,u,v,wrk)
	dimension a(m,n),s(n),u(m,n),v(n,n),wrk(1)
	n2=n+1
	n3=n+n2
	call svd1(a,m,n,m,n,0,n,n,s,u,v,wrk,wrk(n2),wrk(n3))
	return
	end
	subroutine svd1(a,mmax,nmax,m,n,p,nu,nv,s,u,v,b,c,t)
c  singular value decomposition of a matrix
c  Peter A. Businger, Bell Telephone Laboratories
c  Gene H. Golub, Stanford University
c  Algorithm 358, Collected algorithms from ACM
c  modified for real matrix by Mike Webring, USGS
	dimension a(mmax,1), u(mmax,1), v(nmax,1)
	dimension s(n), b(n), c(n), t(n)
	integer p
	data eta,tol/1.5e-8, 1.e-31/
	np=n+p
	n1=n+1
c
c  householder reduction
	c(1)=0.0
	k=1
10	k1=k+1
c
c  elimination of a(i,k), i=k+1...m
	z=0.0
	do 20 i=k,m
20	z = z + a(i,k)**2
	b(k)=0.0
	if(z.le.tol) go to 70
	z=sqrt(z)
	b(k)=z
	w=abs(a(k,k))
	q=1.0
	if(w.ne.0.0) q=a(k,k)/w
	a(k,k)=q*(z+w)
	if(k.eq.np) go to 70
	do 50 j=k1,np
	q=0.0
	do 30 i=k,m
30	q=q+a(i,k)*a(i,j)
	q=q/(z*(z+w))
	do 40 i=k,m
40	a(i,j)=a(i,j)-q*a(i,k)
50	continue
c
c  phase transformation
	q=-sign(1.0,a(k,k))
	do 60 j=k1,np
60	a(k,j)=q*a(k,j)
c
c  elimination of a(k,j), j=k+2...n
70	if(k.eq.n) go to 140
	z=0.0
	do 80 j=k1,n
80	z = z + a(k,j)**2
	c(k1)=0.0
	if(z.le.tol) go to 130
	z=sqrt(z)
	c(k1)=z
	w=abs(a(k,k1))
	q=1.0
	if(w.ne.0.0) q=a(k,k1)/w
	a(k,k1)=q*(z+w)
	do 110 i=k1,m
	q=0.0
	do 90 j=k1,n
90	q = q + a(k,j)*a(i,j)
	q=q/(z*(z+w))
	do 100 j=k1,n
100	a(i,j)=a(i,j)-q*a(k,j)
110	continue
c
c  phase transformation
	q=-sign(1.0,a(k,k1))
	do 120 i=k1,m
120	a(i,k1)=a(i,k1)*q
130	k=k1
	go to 10
c
c  tolerance for negligible elements
140	eps=0.0
	do 150 k=1,n
	s(k)=b(k)
	t(k)=c(k)
150	eps=amax1(eps,s(k)+t(k))
	eps=eps*eta
c
c  initialization of u and v
	if(nu.eq.0) go to 180
	do 170 j=1,nu
	do 160 i=1,m
160	u(i,j)=0.0
170	u(j,j)=1.0
180	if(nv.eq.0) go to 210
	do 200 j=1,nv
	do 190 i=1,n
190	v(i,j)=0.0
200	v(j,j)=1.0
c
c  qr diagonalization
210	do 380 kk=1,n
	k=n1-kk
c
c  test for split
220	do 230 mm=1,k
	m2=k+1-mm
	if(abs(t(m2)).le.eps) go to 290
	if(abs(s(m2-1)).le.eps) go to 240
230	continue
c
c  cancellation of e(m2)
240	cs=0.0
	sn=1.0
	m1=m2-1
	do 280 i=m2,k
	f=sn*t(i)
	t(i)=cs*t(i)
	if(abs(f).le.eps) go to 290
	h=s(i)
	w=sqrt(f*f+h*h)
	s(i)=w
	cs=h/w
	sn=-f/w
	if(nu.eq.0) go to 260
	do 250 j=1,n
	x=u(j,m1)
	y=u(j,i)
	u(j,m1) = x*cs + y*sn
250	u(j,i)  = y*cs - x*sn
260	if(np.eq.n) go to 280
	do 270 j=n1,np
	q=a(m1,j)
	r=a(i,j)
	a(m1,j) = q*cs + r*sn
270	a(i,j)  = r*cs - q*sn
280	continue
c
c  test for convergence
290	w=s(k)
	if(m2.eq.k) go to 360
c
c  origin shift
	x=s(m2)
	y=s(k-1)
	g=t(k-1)
	h=t(k)
	f=((y-w)*(y+w)+(g-h)*(g+h))/(2.0*h*y)
	g=sqrt(f*f+1.0)
	if(f.lt.0.0) g=-g
	f=((x-w)*(x+w)+(y/(f+g)-h)*h)/x
c
c  qr step
	cs=1.0
	sn=1.0
	m1=m2+1
	do 350 i=m1,k
	g=t(i)
	y=s(i)
	h=sn*g
	g=cs*g
	w=sqrt(h*h+f*f)
	t(i-1)=w
	cs=f/w
	sn=h/w
	f=x*cs+g*sn
	g=g*cs-x*sn
	h=y*sn
	y=y*cs
	if(nv.eq.0) go to 310
	do 300 j=1,n
	x=v(j,i-1)
	w=v(j,i)
	v(j,i-1)= x*cs + w*sn
300	v(j,i)  = w*cs - x*sn
310	w=sqrt(h*h+f*f)
	s(i-1)=w
	cs=f/w
	sn=h/w
	f=cs*g+sn*y
	x=cs*y-sn*g
	if(nu.eq.0) go to 330
	do 320 j=1,n
	y=u(j,i-1)
	w=u(j,i)
	u(j,i-1) = y*cs + w*sn
320	u(j,i)   = w*cs - y*sn
330	if(n.eq.np) go to 350
	do 340 j=n1,np
	q=a(i-1,j)
	r=a(i,j)
	a(i-1,j)=q*cs+r*sn
340	a(i,j) = r*cs-q*sn
350	continue
	t(m2)=0.0
	t(k)=f
	s(k)=x
	go to 220
c
c  convergence
360	if(w.ge.0.0) go to 380
	s(k)=-w
	if(nv.eq.0) go to 380
	do 370 j=1,n
370	v(j,k)=-v(j,k)
380	continue
c
c  sort singular values
	do 450 k=1,n
	g=-1.0
	j=k
	do 390 i=k,n
	if(s(i).le.g) go to 390
	g=s(i)
	j=i
390	continue
	if(j.eq.k) go to 450
	s(j)=s(k)
	s(k)=g
	if(nv.eq.0) go to 410
	do 400 i=1,n
	q=v(i,j)
	v(i,j)=v(i,k)
400	v(i,k)=q
410	if(nu.eq.0) go to 430
	do 420 i=1,n
	q=u(i,j)
	u(i,j)=u(i,k)
420	u(i,k)=q
430	if(n.eq.np) go to 450
	do 440 i=n1,np
	q=a(j,i)
	a(j,i)=a(k,i)
440	a(k,i)=q
450	continue
c
c  back transformation
	if(nu.eq.0) go to 510
	do 500 kk=1,n
	k=n1-kk
	if(b(k).eq.0.0) go to 500
	q=-sign(1.0,a(k,k))
	do 460 j=1,nu
460	u(k,j)=q*u(k,j)
	do 490 j=1,nu
	q=0.0
	do 470 i=k,m
470	q = q + a(i,k)*u(i,j)
	q = q/(abs(a(k,k))*b(k))
	do 480 i=k,m
480	u(i,j)=u(i,j)-q*a(i,k)
490	continue
500	continue
510	if(nv.eq.0) go to 570
	if(n.lt.2) go to 570
	do 560 kk=2,n
	k=n1-kk
	k1=k+1
	if(c(k1).eq.0.0) go to 560
	q=-sign(1.0,a(k,k1))
	do 520 j=1,nv
520	v(k1,j)=q*v(k1,j)
	do 550 j=1,nv
	q=0.0
	do 530 i=k1,n
530	q=q+a(k,i)*v(i,j)
	q=q/(abs(a(k,k1))*c(k1))
	do 540 i=k1,n
540	v(i,j)=v(i,j)-q*a(k,i)
550	continue
560	continue
570	return
	end
c******************************************************************************
	subroutine ata(a,m,n,b)
	dimension a(m,n),b(n,n)
	do 10 j=1,n
	do 10 i=1,n
10	b(i,j)=xty(a(1,i),1,a(1,j),1,m)
	return
	end
c******************************************************************************
	subroutine axb(a,ma,na,b,mb,nb,c)
	dimension a(ma,na),b(mb,nb),c(ma,nb)
	if(na.ne.mb) go to 20
	k=na
	do 10 j=1,nb
	do 10 i=1,ma
10	c(i,j)=xty(a(i,1),ma,b(1,j),1,k)
	return
20	print 21
21	format(' axb: col a .ne. row b')
	return
	end
c******************************************************************************
	subroutine atb(a,ma,na,b,mb,nb,c)
	dimension a(ma,na),b(mb,nb),c(na,nb)
	if(ma.ne.mb) go to 20
	k=ma
	do 10 j=1,nb
	do 10 i=1,na
10	c(i,j)=xty(a(1,i),1,b(1,j),1,k)
	return
20	print 21
21	format(' atb: col.ne.row')
	return
	end
c******************************************************************************
	function xty(x,ix,y,iy,n)
	dimension x(n),y(n)
	double precision tmp,tp,tn
	tp=0.0d0
	tn=0.0d0
	jx=1
	jy=1
	do 10 i=1,n
	tmp=dble(x(jx))*dble(y(jy))
	if(tmp.lt.0.0d0) go to 5
	tp=tp+tmp
	go to 6
5	tn=tn+tmp
6	jx=jx+ix
10	jy=jy+iy
	xty=sngl(tp+tn)
	return
	end
c******************************************************************************
	subroutine chki(a,n,iexp)
c  check the identity matrix
	dimension a(n,n)
	iexp=15
	tmx=0.0
	do 10 j=1,n
	do 10 i=1,n
	t=abs(a(i,j))
	if(i.eq.j) t=abs(1.0-a(i,i))
	if(t.gt.tmx) tmx=t
10	continue
	if(tmx.le.1.e-15) return
	iexp=int(0.5-alog10(tmx))
	return
	end
c******************************************************************************
	subroutine fit1(xa,y,n,dx,a,y0)
c  least square fit to  y=ax+y0
	dimension z(5),y(n),xa(1)
	equivalence (z(1),sx),(z(2),sy),(z(3),sxy),
     1 (z(4),sx2),(z(5),x)
	do 1 i=1,5
1	z(i)=0.
	ix=0
	if(dx.eq.0.) ix=1
	a=0.
	y0=y(1)
	if(n.le.1) return
	do 2 i=1,n
	if(ix.gt.0) x=xa(i)-xa(1)
	sx=sx+x
	sy=sy+y(i)
	sxy=sxy+x*y(i)
	sx2=sx2+x*x
2	x=x+dx
	xn=float(n)
	d=xn*sx2-sx*sx
	a=1.e10
	if(abs(d) .lt. 1.e-10) return
	a=(xn*sxy-sx*sy)/d
	y0=(sy-a*sx)/xn
	return
	end
c******************************************************************************
	function correl(n,a,b)
	dimension a(n),b(n)
	correl=0.0
	ab=0.0
	a2=0.0
	b2=0.0
	do 1 i=1,n
	ab=ab+a(i)*b(i)
	a2=a2+a(i)*a(i)
1	b2=b2+b(i)*b(i)
	a2b2=a2*b2
	if(a2b2.lt.1.e-30) return
	correl=abs(ab)/sqrt(a2b2)
	return
	end
c******************************************************************************
	function acosd2(dx,dh)
	acosd2=1.e+37
	hs=dh*dh
	xs=dx*dx
	if(hs.lt.xs) return
	dy=sqrt(hs-xs)
	acosd2=atan2(dy,dx)*57.29578
	return
	end
c******************************************************************************
	logical function inside(n,x,z,xp,zp)
c   is xp,zp inside polygon defined by x,z ?
	dimension x(1),z(1)
	logical z1,z2
	inside=.false.
	do 5 i=1,n
	i1=i+1
	if(i1.gt.n) i1=1
	dz=z(i1)-z(i)
	if(dz.eq.0.0) go to 5
	z1 = zp.le.z(i)
	z2 = zp.gt.z(i1)
	if((z1.and.z2) .or. (.not.z1 .and. .not.z2)) go to 3
	go to 5
c   points on boundary are outside with condition .lt.
c   for inclusion change to .le.
3	rslope=(x(i1)-x(i))/dz
	d=(xp-x(i))-(zp-z(i))*rslope
c        if(abs(d).lt.1.e-4) d=0.0
        if(abs(d).lt.1.e-6) d=0.0
	if(d.lt.0.0) inside = .not.inside
5	continue
	return
	end
c******************************************************************************
	subroutine ammi(n,a,amn,amx,mn,mx,isave)
	dimension a(n)
	if(isave.gt.1) go to 1
	mn=1
	mx=1
	amn=a(1)
	amx=a(1)
1	do 3 i=1,n
	if(a(i).ge.amn) go to 2
	amn=a(i)
	mn=i
2	if(a(i).le.amx) go to 3
	amx=a(i)
	mx=i
3	continue
	return
	end
c******************************************************************************
	logical function member(itest,m,n)
	dimension m(n)
	member=.true.
	do 10 i=1,n
	if(m(i).eq.itest) return
10    continue
	member=.false.
	return
	end
c******************************************************************************
	subroutine shsort(n,w,m,iasend)
c  shell diminishing increment sort
c  iasend < 0 for descending order
	dimension w(n),m(n)
	do 1 i=1,n
1	m(i)=i
	is=n
10	is=is/2
	if(is) 70,70,20
20	k=n-is
	i2=1
30	i=i2
40	j=i+is
	if(iasend) 42,43,43
42	if(w(m(i))-w(m(j))) 60,50,50
43	if(w(m(j))-w(m(i))) 60,50,50
50	i2=i2+1
	if(i2-k) 30,30,10
60	it=m(i)
	m(i)=m(j)
	m(j)=it
	i=i-is
	if(i) 50,50,40
70	return
	end
c******************************************************************************
	subroutine cvc(a,i2)
c  i2=1 for lower case to upper
	dimension ilb(2),iub(2),icv(2)
	character a*(*)
	data ilb/97,65/, iub/122,90/, icv/-32,+32/
	n=len(a)
	if(i2.eq.1 .or. i2.eq.2) then
	do 10 i=1,n
	  j=ichar(a(i:i))
	  if(j.lt.ilb(i2) .or. j.gt.iub(i2)) go to 10
	  j=j+icv(i2)
	  a(i:i)=char(j)
10	continue
	else
	print *,' cvc: invalid convert parameter'
	endif
	return
	end
c******************************************************************************
	function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
	character a*(*)
	n=len(a)
	if(a(1:1).ne.' ') go to 15
	do 1 m=2,n
1	if(a(m:m).ne.' ') go to 5
	leftj=0
	return
c
5	i2=0
	do 10 i=m,n
	i2=i2+1
	a(i2:i2)=a(i:i)
10	continue
	do 11 i3=i2+1,n
11	a(i3:i3)=' '
	n=n-m+1
c
15	do 20 leftj=n,1,-1
20	if(a(leftj:leftj).ne.' ') go to 25
25	return
	end
c******************************************************************************
        function noyes()
	character inp*20
	ic=0
	noyes=-1
1	continue
	read(5,2,err=5) inp
2	format(a20)
	if(ic.ge.3) go to 9
4	if(inp(1:1).eq.'y' .or. inp(1:1).eq.'Y') noyes=1
	if(inp(1:1).eq.'n' .or. inp(1:1).eq.'N') noyes=0
	if(noyes.gt.-1) return
5	print  6
6	format(' y or n:'$)
	ic=ic+1
	go to 1
9	print *,' count exceeded, answering no'
	noyes=0
	return
	end
c******************************************************************************
        subroutine clrscr
c   iop=1 (or ne 2) exit in Tektronix mode,
c   iop=2 exit in ANSI character mode
c   "0g is a Digital Engineering Retrographics command 
c	call sendesc('"0g')
c	call sendesc('[2J')
c	call sendesc('[1;1H')
c	call sendesc(char(29))
c	call sendesc(char(12))
c	if(iop.eq.2) call sendesc('"0g')
c       call setdef
        call setmod(3)
       return
 	end
c******************************************************************************
	subroutine sendesc(tail)
	character fmt*14,esc*1,tail*(*)
	n=len(tail)
	fmt=' '
	assign 1 to label
	if(n.gt.9) assign 2 to label
	if(n.gt.99) assign 3 to label
	write(fmt,label) n
1	format( '(1x,a1,a', i1, ',$)' )	
2	format( '(1x,a1,a', i2, ',$)' )	
3	format( '(1x,a1,a', i3, ',$)' )	
	esc=char(27)
	write(6,fmt) esc,tail
	return
	end
c******************************************************************************
	subroutine symbol(x,y,nxy,ich,size,angle)
c  symbol driver where x,y are in data units (clipping occurs outside
c  the data window).
c  'ich' is an integer akin to fortran 77 coding ich=038=ichar('&').
c  'vchar'  uses 1-13, 24-31 for symbols such as diamonds and squares, 
c  and 32-126 as the normal ascii character sequence when icode=1.
	dimension x(nxy),y(nxy)
	data icode/1/, xoff,yoff/0.0,0.0/
	radian=1.7453292e-2*angle
	call vchar(x,y,ich,nxy,icode,size,radian,xoff,yoff)
	return
	end
c******************************************************************************
	subroutine text(x,y,string,size,angle)
c  text driver where x,y are in plot units.
	dimension is(25)
	character s*100
	character string*(*)
	equivalence (is,s)
	data icode/3/, xoff,yoff/0.0,0.0/
	lenstr=len(string)
	do 5 nch=lenstr,1,-1
	  ich=ichar(string(nch:nch))
	  if(ich.gt.32 .and. ich.lt.127) go to 10
5     continue
	return
10	s(1:nch)=string(1:nch)
	radian=1.7453292e-2*angle
	call vchar(x,y,is,nch,icode+1400,size,radian,xoff,yoff)
	return
	end
c******************************************************************************
      function ran3(ix,iy,iz)
c
c  Following "save" statement commented out because formal arguements
c  can't be saved, but MS Fortran saves all variable by default, anyway.
c  Lin Cordell June 91.
comment out      save ix,iy,iz
c
c     random number generator from algorithm
c     of Brian Wichmann and David Hill
c     in BYTE magazine,March,1987
c
      data ifirst/0/,jx/1/,jy/10000/,jz/3000/
      if(ifirst.eq.0.and.ix.eq.0.and.iy.eq.0
     & .and.iz.eq.0) then
      ifirst=1
      ix=jx
      iy=jy
      iz=jz 
      endif
c
c     first generator
      ixx = ix/177
      ix = 171 * (ix - 177 * ixx) - 2 * ixx
      if(ix.lt.0) ix = ix + 30269
c
c     second generator
      iyy = iy/176
      iy = 172 * (iy - 176 * iyy) - 35 * iyy
      if(iy.lt.0) iy = iy + 30307
c
c     third generator
      izz = iz/178
      iz = 170 * (iz - 178 * izz) - 63 * izz
      if(iz.lt.0) iz = iz + 30323
c
c     combine
      temp = ix/30269. + iy/30307. + iz/30323.
      ran3 = temp - aint(temp)
      return
      end
c******************************************************************************
      subroutine namesk(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,*920)
      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
  920 return 1
      end
c*********
      subroutine check(pvar,tvar,nn,chv,nvar,numa,inum,*)
c
c     assigns values to proper variable
c     variables are passed to program saki through common blocks
c     numr=position in the array var where real variables start
c     numa=position in the array var where arrays start
c     nnvar=number of variables in program minc
c
      parameter(nnvar=50,numr=12)
      character*6 pvar,var(nnvar)
      character*56 tvar,kvar,cfmt
      logical chv
        character compon*1
        character*20 titlx,titlz,titlg,titlh
	character*56 bfile,mfile,gfile,hfile,ifmtg,ifmth,ofmtg,ofmth
	character*56 title1,title2,title3
	common /parm1/  nvert,nbody,nobs,mobs,iswt,jswt,
     1               iave,jave,vex,istatn,iparsh,labelv,numbod
	common /magxyz/ ef(3),ev(3),tv(3),azmuth,icompn
	common /parms/ bfile,mfile,gfile,hfile,ifmtg,ifmth,ofmtg,ofmth,
     1 efield,einc,edec,compon,plenth
	common /plot1/  iplotr,sizea,sizet,sizel,
     1		title1,title2,title3,
     1		xscale,xxx(2),adelx,lintx,titlx,
     1		zscale,zzz(2),adelz,lintz,titlz,
     1		gscale,ggg(2),adelg,lintg,titlg,
     1		hscale,hhh(2),adelh,linth,titlh
	common /pltdev/ jplotr,xboard,yboard,xlimit,ylimit
c
      data var/'iave','jave','istatn','iparsh','labelv','numbod',
     & 'iplotr','lintx','lintz','lintg','linth','efield','einc',
     & 'edec','azmuth','compon','plenth','xlimit','ylimit',
     & 'vex','sizea','sizet','sizel','bfile','mfile','gfile','hfile',
     & 'ifmtg','ifmth','ofmtg','ofmth','title1','title2',
     & 'title3','titlx','titlz','titlg','titlh','xscale',
     & 'zscale','gscale','hscale','adelx','adelz','adelg',
     & 'adelh','xxx','zzz','ggg','hhh'/
      numa=47
      do 190 i=1,nnvar
      if(pvar.ne.var(i)) go to 190
      if(.not.chv) then
c
c     noncharacter value
c     right justify then 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
      go to (101,102,103,104,105,106,107,108,109,110,111,112,
     & 113,114,115,116,117,118,119,120,121,122,123,124,1241,125,
     & 126,127,128,129,130,131,132,133,134,135,136,137,138,
     & 139,140,141,142,143,144,145,146,147,148,149),i
  101 iave=jvar
      go to 200
  102 jave=jvar
      go to 200
  103 istatn=jvar
      go to 200
  104 iparsh=jvar
      go to 200
  105 labelv=jvar
      go to 200
  106 numbod=jvar
      go to 200
  107 iplotr=jvar
      go to 200
  108 lintx=jvar
      go to 200
  109 lintz=jvar
      go to 200
  110 lintg=jvar
      go to 200
  111 linth=jvar
      go to 200
  112 efield=xvar
      go to 200
  113 einc=xvar
      go to 200
  114 edec=xvar
      go to 200
  115 azmuth=xvar
      go to 200
  116 compon=xvar
      go to 200
  117 plenth=xvar
      go to 200
  118 xlimit=xvar
      go to 200
  119 ylimit=xvar
      go to 200
  120 vex=xvar
      go to 200
  121 sizea=xvar
      go to 200
  122 sizet=xvar
      go to 200
  123 sizel=xvar
      go to 200
  124 bfile=tvar(1:nn)
      go to 200
 1241 mfile=tvar(1:nn)
      go to 200
  125 gfile=tvar(1:nn)
      go to 200
  126 hfile=tvar(1:nn)
      go to 200
  127 ifmtg=tvar(1:nn)
      go to 200
  128 ifmth=tvar(1:nn)
      go to 200
  129 ofmtg=tvar(1:nn)
      go to 200
  130 ofmth=tvar(1:nn)
      go to 200
  131 title1=tvar(1:nn)
      go to 200
  132 title2=tvar(1:nn)
      go to 200
  133 title3=tvar(1:nn)
      go to 200
  134 titlx=tvar(1:nn)
      go to 200
  135 titlz=tvar(1:nn)
      go to 200
  136 titlg=tvar(1:nn)
      go to 200
  137 titlh=tvar(1:nn)
      go to 200
  138 xscale=xvar
      go to 200
  139 zscale=xvar
      go to 200
  140 gscale=xvar
      go to 200
  141 hscale=xvar
      go to 200
  142 adelx=xvar
      go to 200
  143 adelz=xvar
      go to 200
  144 adelg=xvar
      go to 200
  145 adelh=xvar
      go to 200
  146 xxx(inum)=xvar
      go to 200
  147 zzz(inum)=xvar
      go to 200
  148 ggg(inum)=xvar
      go to 200
  149 hhh(inum)=xvar
      go to 200
  190 continue
      write(*,*) ' error in namelist - ',pvar,' variable not included'
      return 1
  200 nvar=i
      return
      end
 
