c                  USGS Open File-Report 85-122
c
c                      program saki
c
c    Semi-Automatic MarQuardt Inversion of gravity and magnetic 
c         profiles using singular value decomposition.
c
c  computer: DEC VAX
c  language: fortran 77
c  development history:
c  original gravity inversion program 10/80
c  2.5d gravity and magnetics 11/81
c  simultaneous inversion 6/83
c
c      Mike Webring
c  Br. of Geophysics, mail stop 964
c  US Geological Survey, box 25046
c  Denver Federal Center
c  Lakewood, CO  80225
c
        common /field/  nf(10),f(5120)
	common /model/  p(900)
	common /vertx/  iv(1001)
	common /swork/  a(8420)
	common /swork2/  a2(8400)
	common /parms/ bfile,mfile,gfile,hfile,ifmtg,ifmth,ofmtg,ofmth,
     1 efield,einc,edec,compon,plenth
	common /parm1/  nvert,nbody,nobs,mobs,iswt,jswt,
     1               iave,jave,vex,istatn,iparsh,labelv,numbod
	common /parm2/  iofmt(14,4)
	common /magxyz/ ef(3),ev(3),tv(3),azmuth,icompn
	common /freel/  weight,rms1,iter,xk1,nfree(4),ifree(80)
	common /btags/  itag(450)
	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
	common /window/ xwind(2,4),ywind(2,4)
	common /vuport/ xview(4,4),yview(4,4)
	common /scal1/  onport(12)
c
      dimension iaa(8420)
	dimension xb(200),zb(200),rho(50),sus(50),pl1(50),pl2(50),
     1 rems(50),remi(50),remd(50)
	dimension icrn(500),ncrn(150),loca(150),ibord(200)
	dimension ipar(9),iparm(9)
c
	character    bfile*56, ifun*20,compon*1,iparm*4,prog*8
        character*20 titlx,titlz,titlg,titlh
	character*56 test,cf,mfile,gfile,hfile,ifmtg,ifmth,ofmtg,ofmth
	character*56 title1,title2,title3
       character*56 jfmtg,jfmth,kfmtg,kfmth
c
      equivalence(a(1),iaa(1))
	equivalence (iswt,grvswt), (jswt,magswt)
	equivalence (p(1),xb), (p(201),zb), (p(401),rho), (p(451),sus),
     1 (p(501),pl1), (p(551),pl2), (p(601),rems), (p(651),remi),
     1 (p(701),remd)
	equivalence (iv(1),icrn),    (iv(501),ncrn), (iv(651),loca),
     1           (iv(801),nbord), (iv(802),ibord)
       equivalence (iofmt(1,1),jfmtg), (iofmt(1,2),jfmth),
     1           (iofmt(1,3),kfmtg), (iofmt(1,4),kfmth)

c
	logical grvswt,magswt
c
	data npar/9/,ipar/1,201,401,451,501,551,601,651,701/
	data iparm/'xbod','zbod','rho ','susc','pl1 ','pl2 ','rems',
     1 'remi','remd'/
        data maxfld/5120/
c
	prog='saki'
	idum=0
	ioutp=0
	test=' '
	mfile=' '
	bfile=' '
	gfile=' '
	hfile=' '
	ifmtg=' '
	ifmth=ifmtg
	ofmtg='(1p3g15.5)'
	ofmth=ofmtg
	compon=' '
	icompn=0
	efield=50000.0
	einc=90.
	edec=0.
	iave=1
	jave=1
	vex=1.0
	xxx(1)=0.0
	xxx(2)=0.0
	istatn=1
	iparsh=0
	labelv=0
	numbod=1
        do 1 i=1,maxfld
1	f(i)=0.
	do 2 i=1,900
2	p(i)=0.
	nxf=0
	nzf=0
	nrf=0
	call initp
c
109	print  110
110	format(' enter command filename :'$)
	read(5,37,end=109) cf
	if(cf.ne.test) then
	  open (unit=9,file=cf,status='old',form='formatted')
	  call namesk(9,*111) 
	  if(bfile(1:1).ne.' ') mfile=bfile
          jfmtg=ifmtg
          jfmth=ifmth
          kfmtg=ofmtg
          kfmth=ofmth
          ef(1)=efield
          ef(2)=einc
          ef(3)=edec
	  jplotr=iplotr
  	  go to 112
c111	  print *,' there is either a syntax error in the command file'
c	  print *,' or and invalid parameter.'
111	  print *,' do you want a list of valid namelist parameters ?'
	  if(noyes(ldum).eq.1) then
	    print *,' mfile,gfile,hfile,ifmtg,ifmth,ofmtg,ofmth'
	    print *,' iave,jave,efield,einc,edec,azmuth,compon,plenth'
	    print *,' istatn,iparsh,labelv,numbod'
	    print *,' iplotr,xlimit,ylimit,vex,sizea,sizet,sizel'
	    print *,' title1*56,title2*56,title3*56'
	    print *,' xscale,xxx(2),adelx,lintx,titlx*20'
	    print *,' zscale,zzz(2),adelz,lintz,titlz*20'
	    print *,' gscale,ggg(2),adelg,lintg,titlg*20'
	    print *,' hscale,hhh(2),adelh,linth,titlh*20'
	    print *,' note bfile has been changed to mfile'
	  endif
 	  close(9)
 	  stop
112	  close(9)
	endif
	if(mfile.eq.test) then
	  write(*,114)
114      format(' model filename: '$)
	  read(5,37) mfile
	endif
	if( gfile.eq.test) then
	  write(*,115)
115      format(' gravity data filename (optional): '$)
	  read(5,37) gfile
	endif
	if(hfile.eq.test) then
	  write(*,116)
116      format(' magnetic data filename (optional): '$)
	  read(5,37) hfile
	endif
c
c  open data files, set function switches
	grvswt=.false.
	magswt=.false.
	open(unit=10,file=mfile,status='old',form='formatted')
	if(gfile.eq.test) go to 120
	grvswt=.true.
	open(unit=11,file=gfile,status='old',form='formatted',
     1 blank='zero')
120	if(hfile.eq.test) go to 121
	magswt=.true.
	open(unit=12,file=hfile,status='old',form='formatted',
     1 blank='zero')
c
c  read data
121	nobs=0
	if(.not.grvswt) go to 131
125	read(11,*,end=130)
	nobs=nobs+1
	go to 125
130	rewind 11
131	mobs=0
	if(.not.magswt) go to 141
135	read(12,*,end=140)
	mobs=mobs+1
	go to 135
140	rewind 12
c  place holders to avoid dimension errors
141	if(nobs.eq.0) nobs=1
	if(mobs.eq.0) mobs=1
c  data is stored in field common :xg,zg,go,gc,gerr, xh,zh,ho,hc,herr
	call initnf(nobs,mobs,nf,maxfld,ierr)
	if(ierr.eq.1) stop
	call datinp(nobs,f(nf(1)),f(nf(2)),f(nf(3)),
     1	    mobs,f(nf(6)),f(nf(7)),f(nf(8)))
c
c  final initialization
	if (xxx(1).ge.xxx(2)) call initw(xxx,grvswt,magswt)
	call initl(plenth,pl,nbody,0)
	call cvc(compon,2)
	if(compon.eq.'x') icompn=1
	if(compon.eq.'y') icompn=2
	if(compon.eq.'z') icompn=3
	if(magswt) call initm(0)
	call outerb(nvert,nbody,iaa,iaa(301),iaa(601))
	call ckmodl(2,1,dum1,dum2,ierr1)
	if(grvswt) call ckmodl(1,nobs,f(nf(1)),f(nf(2)),ierr2)
	if(magswt) call ckmodl(1,mobs,f(nf(6)),f(nf(7)),ierr2)
	if(.not.grvswt .and. .not.magswt) go to 290
c
c  function branches
c
200	print  201
201	format(' function:'$)
	read(5,37,end=200) ifun
37	format(a20)
	call cvc(ifun,2)
	if(ifun(1:2).eq.'ca') go to 205
	if(ifun(1:2).eq.'in') go to 210
	if(ifun(1:2).eq.'wi') go to 220
	if(ifun(1:2).eq.'ty') go to 230
	if(ifun(1:2).eq.'pl') go to 240
	if(ifun(1:2).eq.'ed') go to 250
	if(ifun(1:2).eq.'ou') go to 260
	if(ifun(1:2).eq.'is') go to 270
	if(ifun(1:2).eq.'se') go to 280
	if(ifun(1:2).eq.'gr') go to 295
	if(ifun(1:2).eq.'ma') go to 295
	if(ifun(1:2).eq.'bo') go to 295
	if(ifun(1:2).eq.'sy') go to 300
	if(ifun(1:2).eq.'co') go to 310
	if(ifun(1:2).eq.'tr') go to 330
	if(ifun(1:2).eq.'no') go to 340
	if(ifun(1:2).eq.'ex') go to 999
	if(ifun(1:2).eq.'st') go to 999
	print  103
103	format(' help ? '$)
	if(noyes(l).eq.0) go to 200
	print  102
102	format(' functions: calculate, invert, plot, edit, print ',/,
     1 ' setup, window, synthetic, isolate body, trend removal ',/,
     1 ' noise, gravity, magnetic, both, output, exit')
	go to 200
c
c  forward calculation
205	call taldvr(nobs,f(nf(1)),f(nf(2)),f(nf(3)),f(nf(4)),f(nf(5)),
     1 rmsg,mobs,f(nf(6)),f(nf(7)),f(nf(8)),f(nf(9)),f(nf(10)),rmsh)
	if(grvswt) print 206,rmsg
206	format('  gravity rms error ',1pg13.3)
	if(magswt) print 207,rmsh
207	format(' magnetic rms error ',1pg13.3)
	go to 200
c
c  inversion
210	call nlctl(nobs,f(nf(1)),mobs,f(nf(6)))
	ngrv=0
	nmag=0
	if(grvswt) ngrv=nobs
	if(magswt) nmag=mobs
	no=ngrv+nmag
	go to 200
c
c  change operating window
220	print *,' enter profile min, max (0,0 to default) :'
	read(5,*,end=200,err=200) xxx
	if(xxx(1).ge.xxx(2)) call initw(xxx,grvswt,magswt)
	go to 200
c
c  print model parameter values
230	call getadr(npar,ipar,iparm,iad)
	call typdat(p(iad))
	go to 200
c
c  plot data and model
240	call pltdvr(ier)
	go to 200
c
c  edit model parameters
250	call getadr(npar,ipar,iparm,iad)
	call edtdat(p(iad))
 	go to 200
c
c  output model and/or data
260	call modsav(ofmtg,ofmth)
	ioutp=1
	go to 200
c
c  plot isolated body curves
270	mgswt=0
	if(magswt) mgswt=1
	if(grvswt.and.magswt) mgswt=-1
	call setwin(mgswt,xxx,zzz,ggg,hhh)
	call isodvr
	go to 200
c
c  change parameters
280	call setup
	go to 200
c
c  switch calculations
290	print 291
291	format(' activate gravity, magnetic, or both :')
	read(5,37,end=290) ifun
	call cvc(ifun,2)
295	grvswt=.false.
	magswt=.false.
	if(ifun(1:2).eq.'gr') grvswt=.true.
	if(ifun(1:2).eq.'ma') magswt=.true.
	if(ifun(1:2).eq.'bo') grvswt=.true.
	if(ifun(1:2).eq.'bo') magswt=.true.
	iter=1
	if(.not.grvswt .and. .not.magswt) go to 290
	if(magswt) call initm(0)
	if(grvswt .and. nobs.le.1) go to 300
	if(magswt .and. mobs.le.1) go to 300
	go to 200
c
c  synthetic data
c  field values are stored xg,zg,go,gc,gerr, xh,zh,ho,hc,herr
300	if(.not.grvswt) go to 305
	if(nobs.gt.1) then
	  print 411
411	  format(' change gravity locations ?'$)
	  if(noyes(ldum).eq.0) go to 305
	endif
	print  301
301	format( ' enter synthetic gravity   xo,dx,nobs,z :'$)
	read(5,*,end=305,err=305) x,dx,nobs,zl
	if(nobs.le.0) nobs=1
c  move magnetic data
	 iad=nf(6)
	 do 400 i=1,5*mobs
	 a(i)=f(iad)
400	 iad=iad+1
	 do 401 i=1,5*nobs
401	 f(i)=0.0
	 iad=5*nobs+1
	 do 402 i=1,5*mobs
	 f(iad)=a(i)
402	 iad=iad+1
	call initnf(nobs,mobs,nf,maxfld,ierr)
	if(ierr.eq.1) go to 300
	print *,'is the z level a constant clearance distance ?'
	if (noyes(ldum).eq.0) then
	  do 403  i = 0, nobs-1
	    f(nf(1)+i) = x
	    f(nf(2)+i) = zl
	    x = x + dx
403       continue	  
	else
	  call drpsur(x,dx,nobs,zl, f(nf(1)),f(nf(2)) )
	endif
c  enter synthethic mag locations
305	if(.not.magswt) go to 309
	if(mobs.gt.1) then
	  print 406
406	  format(' change magnetic locations ?'$)
	  if(noyes(l)) 309,309,410
	endif
410	print  306
306	format( ' enter synthetic magnetic  xo,dx,nobs,z :'$)
	read(5,*,end=309,err=309) x,dx,mobs,zl
	if(mobs.le.0) mobs=1
	call initnf(nobs,mobs,nf,maxfld,ierr)
	if(ierr.eq.1) go to 410
	print *,'is the z level a constant clearance distance ?'
	if (noyes(ldum).eq.0) then
	  do 307 i = 0, mobs-1
	    f(nf(6)+i) = x
	    f(nf(7)+i) = zl
	    x = x + dx
307       continue	   
	else
	  call drpsur(x,dx,mobs,zl, f(nf(6)),f(nf(7)) )
	endif
309	continue
	call taldvr(nobs,f(nf(1)),f(nf(2)),f(nf(3)),f(nf(4)),f(nf(5)),
     1 rmsg,mobs,f(nf(6)),f(nf(7)),f(nf(8)),f(nf(9)),f(nf(10)),rmsh)
c
c  copying function
310	if(.not.grvswt .or. nobs.le.1) go to 314
	print  311
311	format(' copy calculated gravity into observed ?'$)
	if(noyes(l).eq.0) go to 314
	do 313 i=0,nobs-1
313	f(nf(3)+i)=f(nf(4)+i)
c
314	if(.not.magswt .or. mobs.le.1) go to 200
	print  315
315	format(' copy calculated magnetics into observed ?'$)
	if(noyes(l).eq.0) go to 200
	do 316 i=0,mobs-1
316	f(nf(8)+i)=f(nf(9)+i)
	go to 200
c
c  trend removal
330	if(grvswt) then
	call trend(f(nf(1)),f(nf(3)),nobs,dc,dy)
	call ammi(nobs,f(nf(3)),gmin,gmax,mn,mx,1)
	print  335,gmin,gmax
335	format(' residual data is stored in the observation array.',/,
     1 ' the new gravity data has  a min/max of :',1p2e15.5)
	endif
	if(magswt) then
	call trend(f(nf(6)),f(nf(8)),mobs,dc,dy)
	call ammi(mobs,f(nf(8)),gmin,gmax,mn,mx,1)
	print  336,gmin,gmax
336	format(' residual data is stored in the observation array.',/,
     1 ' the new magnetic data has  a min/max of :',1p2e15.5)
	endif
	go to 200
c
c  add noise to data
340	call adnois(nf,f,a,grvswt,magswt)
	go to 200
c
999	close(10)
	if(ioutp.eq.0) then
	  print *,' write the current model to disk ?'
	  if(noyes(ldum).eq.1) call modsav(ofmtg,ofmth)
	endif
	stop
	end
c******************************************************************************
	subroutine datinp(nobs,xg,zg,go, mobs,xh,zh,ho)
c  example model contains vertex coordinate pairs followed by
c  2 line body specifications, the first line is physical parameters
c  and the second line a lookup list which specifies vertex number in
c  clockwise order around the body.
c  -2 1  2 1  2 2  -2 2
c  0 1 0 3  0 -1
c  <<<<
c  .2 1.e-3 20. -20. 0.e-5 40. 0.
c  1 5 6 4 <
c  -.2 -1.e-3 20. -20. 0.e-5 60. 0.
c  5 2 3 6 <
	common /model/ xb(200),zb(200),rho(50),sus(50),
     1 pl(50,2),rem(50,3),dum(150)
	common /vertx/ icrn(500),ncrn(150),loca(150),dum1(201)
	common /parm1/ nvert,nbody,nobs1,mobs1,grvswt,magswt,
     1 iave,jave,vex,idepth,iparsh,indiv,idum
	common /parm2/ fmtg,fmtm,dum3(28)
 	common /swork/ tx(400),tz(400),tg(400),link(400),dum4(6820)
	common /magxyz/ ef(3),ev(3),tv(3),azmuth,idum1
	dimension xg(1),zg(1),go(1),xh(1),zh(1),ho(1)
	logical grvswt,magswt,nofmt
	character*56 fmtg,fmtm
	character*80 line,blank
 	blank=' '
c
	iunit=10
	call countv(nvert,iunit,iend)
	if(iend.eq.1) stop ' datinp: eof1 model file'
	nbody=0
4	nbody=nbody+1
5	read(iunit,'(a)',end=6) line
       if(line.eq.blank) go to 5
	call countv(ncrn(nbody),iunit,iend)
	if(iend.eq.1) go to 6
	go to 4
6	nbody=nbody-1
	rewind 10
	nvert=nvert/2
	print 7,nvert,nbody
7	format(' number of vertices =',i4,/,
     1 ' number of bodies   =',i4)
c  read vertex coordinates
	read(10,*,err=66) (xb(i),zb(i),i=1,nvert)
c  skip blank lines, find body list delimiter
	do 10  i=1,5
	  read(10,*,err=11) delimt
10      continue	 
c  read physical property and vertex lookup lists
c  put in read for pc Microsoft Fortran
11     read(10,'(a)') line
  	j=1
	ivert=0
	is=1
	do 20 ibody=1,nbody
	loca(ibody)=is-1
	ie=loca(ibody)+ncrn(ibody)
c  blank line before next body spec?
	read(10,13) line
13	format(a)
	if (line.ne.blank) backspace (10)
	read(10,*,err=77,end=99) rho(ibody),sus(ibody),pl(ibody,1),
     1 pl(ibody,2), rem(ibody,1),rem(ibody,2),rem(ibody,3)
	read(10,*,err=88,end=99) (icrn(i),i=is,ie)
	is=ie+1
20	continue
c
c  get body tags
25	call bodyid(nvert,nbody)
c
c  read field observations
	if(.not.grvswt) go to 40
	nofmt=.true.
	if( fmtg(1:1).eq.'(' ) nofmt=.false.
	do 30 i=1,nobs
	if(nofmt) then
	read(11,*) xg(i),zg(i),go(i)
	else
	read(11,fmtg) xg(i),zg(i),go(i)
	endif
30	continue
40	if(.not.magswt) go to 160
	nofmt=.true.
	if( fmtm(1:1).eq.'(' ) nofmt=.false.
	do 50 i=1,mobs
	if(nofmt) then
	read(12,*) xh(i),zh(i),ho(i)
	else
	read(12,fmtm) xh(i),zh(i),ho(i)
	endif
50	continue
c
c  sort gravity observations in x
160	if(nobs.eq.0 .and. mobs.eq.0) return
	do 161 i=2,nobs
	if(xg(i).lt.xg(i-1)) go to 165
161	continue
	return
165	call shsort(nobs,xg,link,0)
	do 166 i=1,nobs
	tx(i)=xg(i)
	tz(i)=zg(i)
166	tg(i)=go(i)
	do 167 i=1,nobs
	m=link(i)
	xg(i)=tx(m)
	zg(i)=tz(m)
167	go(i)=tg(m)
	return
66	print *,' datinp: vertex read error'
	return
77	print *,' datinp: phys-prop input error body number',ibody
	return
88	print *,' datinp: look-up read error body number',ibody
	return
99	print *,' datinp: eof in body definition'
	return
	end
c******************************************************************************
	subroutine countv(ival,iunit,iend)
c  count words delimited by spaces or commas
c  return when non-alpanumeric character begins word
	character a*200
	iend=0
	ival=0
10	a=' '
	read(iunit,20,end=99) a
20	format(a200)
30	na=leftj(a)
	if(na.eq.0) go to 10
	if(a(1:1).eq.',') then
	a(1:na-1)=a(2:na)
	a(na:na)=' '
	go to 30
	endif
	ich=ichar(a(1:1))
	if(ich.le.42 .or. ich.ge.123) return
	if(ich.ge.58 .and. ich.le.64) return
	if(ich.ge.91 .and. ich.le.96) return
	mb=index(a,' ')
	mc=index(a,',')
	m=mb
	if(mc.gt.0 .and. mc.lt.mb) m=mc	
	ival=ival+1
	na1=na-m+1
	a(1:na1)=a(m:na)
	do 40 i=na1+1,200
40	a(i:i)=' '
	go to 30
99	iend=1
	return
	end
c******************************************************************************
	subroutine bodyid(nvert,nbody)
c  get the laundry tag at the end of the body lookup list
c        common /swork/ a(2000),dum(6420)
        common /swork/ a(5120),dum(3300)
	common /btags/ itag(450)
	character*36 tag(50)
	character line*80,blank*80,ischar*1
	equivalence (itag,tag)
	logical alphan
	data blank/' '/
	lenl=80
	lentag=36
	rewind 10
	read(10,*,err=5,end=99) (a(i),i=1,2*nvert)
	do 3 i=1,5
	  read(10,*,err=5,end=99) delimt
3       continue
c  put in read for pc Microsoft Fortran
5      read(10,'(a)') line
 	do 100 ibody=1,nbody
	tag(ibody)=' '
10	read(10,20,err=99,end=99) line
	if (line.eq.blank) go to 10
15  	line=' '
	read(10,20,err=99,end=99) line
20	format(a80)
	do 50 ich=1,lenl
	if(.not.alphan(line(ich:ich))) then
	  do 40 is=ich+1,lenl
	  ischar=line(is:is)
	  if(.not.(alphan(ischar)) .or. ischar.eq.' ') go to 40
	  ie = is + (lentag-1)
	  if(ie.gt.lenl) ie=lenl
	  tag(ibody)=line(is:ie)
	  go to 100
40	  continue
	endif
50	continue
	go to 15
100	continue
99	return
	end
c******************************************************************************
	logical function alphan(ch)
c  is 'ch' alpha-numeric ?
	character ch*(1)
	alphan=.false.
	i=ichar(ch)
	if(i.eq.32 .or.  i.eq.43 .or. i.eq.44)  go to 10
	if(i.eq.45 .or.  i.eq.46)  go to 10
	if(i.ge.48 .and. i.le.57)  go to 10
	if(i.ge.65 .and. i.le.90)  go to 10
	if(i.ge.97 .and. i.le.122) go to 10
	return
10	alphan=.true.
	return
	end
c******************************************************************************
	subroutine outerb(nvert,nbody,idupl,is,ie)
c   assemble an array defining the outer boundary
c   of the model for checking vertex relations.
	common /vertx/ icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	dimension idupl(1),is(1),ie(1)
	data maxs/300/,maxb/200/
	nbord=0
	ic=1
	do 10 j=1,nbody
	do 10 i=1,ncrn(j)
	idupl(ic)=0
	is(ic)=icrn(loca(j)+i)
	i1=i+1
	if(i.eq.ncrn(j)) i1=1
	ie(ic)=icrn(loca(j)+i1)
	if((is(ic).lt.1 .or. is(ic).gt.nvert) .or.
     1  (ie(ic).lt.1 .or. ie(ic).gt.nvert)) then
	  print *,' outerb:  lookup list for body ',j
	  print *,' contains an entry < zero or >', nvert, '(nvert)'
	  print *,' corners',i,i1,' point to vertices',is(ic),ie(ic)
	  stop
	endif
	ic=ic+1
	if(ic.gt.maxs) go to 999
10	continue
	nseg=ic-1
c
	do 50 iseg=1,nseg
	iv1=is(iseg)
	iv2=ie(iseg)
	do 20 itst=iseg+1,nseg
	if(iv1.eq.is(itst) .and. iv2.eq.ie(itst)) go to 15
	if(iv1.eq.ie(itst) .and. iv2.eq.is(itst)) go to 15
	go to 20
15	idupl(iseg)=idupl(iseg)+1
	idupl(itst)=idupl(itst)+1
20	continue
50	continue
c
	icount=0
	do 60 i=1,nseg
	if(idupl(i).gt.1) then
	print  61,is(i),ie(i)
61	format(' segment', 2i4,' found more than twice')
	icount=icount+1
	if(icount.gt.5) stop
	endif
	if(idupl(i).eq.0) idupl(i)=-1
60	continue
c
	ipass=0
	do 70 ifirst=1,nseg
	if(idupl(ifirst).ge.0) go to 70
	ibord(1)=is(ifirst)
	ibord(2)=ie(ifirst)
	nbord=2
	go to 71
70	continue
c
71	ifind=ifirst+1
	ipass=ipass+1
	if(ipass.gt.nseg) go to 777
72	if(idupl(ifind).ge.0) go to 75
	if(is(ifind).ne.ibord(nbord)) go to 75
	nbord=nbord+1
	if(nbord.gt.maxb) go to 999
	ibord(nbord)=ie(ifind)
	if(ibord(nbord).eq.ibord(1)) go to 80
	idupl(ifind)=0
	go to 71
75	ifind=ifind+1
	if(ifind.gt.nseg) go to 71
	go to 72
80	nbord=nbord-1
	return
777	print *,' outerb: did not close polygon'
	print *,' current vertex array'
	print *,(ibord(i),i=1,nbord)
999	print  998
998	format(' no model boundary check will be done')
	return
	end
c******************************************************************************
	subroutine ckmodl(ifun,npt,xf,zf,ierr)
c  ifun=0 check both terrain and prisms
c  ifun=1 check terrain, ifun=2 check prism twists
	common /parm1/ nvert,nbody,nobs,mobs,grvswt,magswt,dum(7)
	common /model/ xb(200),zb(200),dum1(500)
	common /vertx/icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	common /swork/ xf1(400),zf1(400),x(100),z(100),dum2(7420)
	dimension lint(2),xf(npt),zf(npt)
	logical grvswt,magswt,inside,member,planar
	ierr=0
	nn=loca(nbody)+ncrn(nbody)
	if(ifun.eq.2) go to 35
c
c   check for vertices above terrain (observations)
	do 20 i=1,npt
	xf1(i)=xf(i)
20	zf1(i)=zf(i)
	m=npt+2
	xf1(npt+1)=xf1(npt)
	xf1(m)=xf1(1)
	zf1(npt+1)=-100.
	zf1(m)=-100.
	do 30 i=1,nvert
	if(.not.member(i,icrn,nn)) go to 30
	if(inside(m,xf1,zf1,xb(i),zb(i))) print  25,i
25	format(' vertex',i3,' is above observations')
30	continue
	if(ifun.eq.1) return
c
35	do 50 j=1,nbody
	if(planar(ncrn(j),xb,zb,icrn(loca(j)+1),lint)) go to 50
	ierr=1
	print *,' ckmodl>> body number ',j,' has an intersection'
	print *,'       between line segments', lint(1), 'and', lint(2)
50	continue
	if(planar(nbord,xb,zb,ibord,lint)) return
	ierr=1
	print *,' ckmodl>> the main model outer boundary has an'
	print *,'          intersection between line segments'
	print *,lint(1), 'and', lint(2)
	return
	end
c******************************************************************************
	subroutine getadr(npar,ipar,iparm,iad)
	dimension ipar(npar),iparm(npar)
	character what*8,iparm*4
30	print  31
31	format(' specify parameter :'$)
	read(5,5,end=30) what
5	format(a8)
	call cvc(what,2)
	do 32 i=1,npar
	if(what(1:4).eq.iparm(i)) go to 34
32	continue
	print  33,iparm
33	format(' model parameters available (enter 4 char.)',/,
     1 10(1x,a4,',') )
	go to 30
34	iad=ipar(i)
	return
	end
c******************************************************************************
	subroutine typdat(a)
	dimension a(1)
1	print  2
2	format(' start,end :'$)
	read(5,*,end=20,err=20) is,ie
	if(is.le.0) return
	if(ie-is+1 .gt. 50) then
	  print *,' limiting output to 50 values'
	  ie=is+49
	endif 
	print  10,(a(i),i=is,ie)
10	format(5g13.5)
	go to 1
20	return
	end
c******************************************************************************
	subroutine edtdat(a)
	dimension a(1)
	character func*16
	ierr=0
1	print  2
2	format(' ed, ty, ret : '$)
	read(5,3,end=1) func
3	format(a16)
	call cvc(func,2)
	nch=leftj(func)
	if(func(1:1).ne.'e') go to 10
	  print 5	
5	  format(' enter n pairs (index1, value1)...'$)
	  read(5,*,end=1,err=1) n,(i,a(i),m=1,n)
	  go to 1
10	if(func(1:1).ne.'t') go to 15
	  call typdat(a)
	  go to 1
15	if(func(1:1).eq.'r') return
	  ierr=ierr+1
	  if(ierr.gt.3) return
	  print 20
20	  format(' edit functions = edit, print, return')
	go to 1
	end
c******************************************************************************
	subroutine modsav(ofmtg,ofmtm)
        common /field/  nf(10),f(5120)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,dum(7)
	common /model/  xb(200),zb(200),rho(50),sus(50),
     1               pl(50,2),rem(50,3),dum1(150)
	common /vertx/  icrn(500),ncrn(150),loca(150),idum(201)
	common /swork/  a(8420)
	common /btags/  itag(450)
	logical grvswt,magswt
        character ofmtg*56,ofmtm*56,ians*1,fmt1*24
	character*36 tag(50),blk20
	equivalence (itag,tag)
	blk20=' '
37      format(a)
c
	open(unit=12,file='model.out',form='formatted',
     1 status='unknown')
	ic=0
	do 2 i=1,nvert
	ic=ic+1
	a(ic)=xb(i)
	ic=ic+1
2	a(ic)=zb(i)
	n2=2*nvert
	write(12,5,err=8) (a(i),i=1,n2)
c  5 lines (20 vertices) then a blank line
5       format( 5( 5(1x, 4(f8.3,1x,f8.3,2x),/) ,/) )
	go to 9
8	print *,' modsav:lost some vertices due to formatting error'
c
9	write(12,10,err=15)
10	format(' >> body specifications')
c
15	do 50 i=1,nbody
	write(12,20,err=21) rho(i),sus(i),pl(i,1),pl(i,2),
     1 rem(i,1),rem(i,2),rem(i,3),i
20	format((1x,f6.3), (1x,1pg11.3), 0p2(1x,f8.1),
     1 (1x,1pg11.3), 0p2(1x,f6.1), ' <<phys_prop',i2)
c
21	nleft=ncrn(i)
	if(nleft.gt.200) go to 50
	iadr=loca(i)
	is=1
25	if(nleft.le.10) go to 35
	write(12,30) (icrn(i2),i2=iadr+is,iadr+is+9)
30	format(10i4)
	nleft=nleft-10
	is=is+10
	go to 25
c
35	write(fmt1,40,err=45) nleft
40	format( '(', i2, 'i4,'' < '',a36)')
	if(tag(i)(1:36).eq.blk20(1:36)) then
	  write(tag(i),43) i
43	  format('body',i3)
	endif
	write(12,fmt1,err=45) 
     1 (icrn(i2),i2=iadr+is,iadr+ncrn(i)),tag(i)
	go to 50
45	print *,' modsav:lost a lookup list for body',i 	
50	continue
	close(12)
c
c  output field values
	print 51
51	format(' output field values ?')
	if(noyes(ldum).eq.0) return
55	print 60
60	format(' calculated or observed :')
	read(5,37,end=55) ians
	call cvc(ians,2)
	ical=-1
	if(ians(1:1).eq.'o') ical=0
	if(ians(1:1).eq.'c') ical=1
	if(ical.eq.-1) go to 55
c
	if(grvswt) then
	open(unit=13,file='gfield.out',form='formatted',
     1 status='unknown')
	n1=nf(1)
	n2=nf(2)
	n3=nf(3+ical)
	do 70 i=0,nobs-1
70	write(13,ofmtg,err=71) f(n1+i),f(n2+i),f(n3+i)
71	close(13)
	endif
c
	if(magswt) then
	open(unit=13,file='hfield.out',form='formatted',
     1 status='unknown')
	n1=nf(6)
	n2=nf(7)
	n3=nf(8+ical)
	do 90 i=0,mobs-1
	write(13,ofmtm,err=91) f(n1+i),f(n2+i),f(n3+i)
90	continue
91	close(13)
	endif
	return
	end
c******************************************************************************
	subroutine taldvr(nobs,xg,zg,go,gc,gerr,rmsg,
     1 mobs,xh,zh,ho,hc,herr,rmsm)
c   basic forward calculation, the talwani driver cycles through
c   the body array, removes a dc level, and figures rms error.
	common /parm1/ nvert,nbody,nobs1,mobs1,grvswt,magswt,
     1 iave,jave,vex,idepth,iparsh,indiv,idum
	common /swork/ xb1(200),zb1(200),dum(8020)
	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 xg(1),zg(1),go(1),gc(1),gerr(1),
     1 xh(1),zh(1),ho(1),hc(1),herr(1)
	logical inside,grvswt,magswt
c
	if(grvswt) then
	do 1 i=1,nobs
1	gc(i)=0.0
	rmsg=1.e30
	endif
	if(magswt) then
	do 2 i=1,mobs
2	hc(i)=0.0
	rmsm=1.e30
	endif
	call ckmodl(2,nobs,xg,zg,ierr)
c	if(grvswt) call ckmodl(1,nobs,xg,zg,ierr)
	if(magswt) call ckmodl(1,mobs,xh,zh,ierr)
	if(ierr.eq.1) return
c
c  check for overlapping bodies
	ic=1
	do 100 j=1,nbody
	do 5 i=1,ncrn(j)
	xb1(i)=xb(icrn(loca(j)+i))
5	zb1(i)=zb(icrn(loca(j)+i))
	do 6 k=1,nvert
	do 9 l=1,ncrn(j)
	if(k.eq.icrn(loca(j)+l)) go to 6
9	continue
	if(.not.inside(ncrn(j),xb1,zb1,xb(k),zb(k))) go to 6
	print  7,k,j
7	format(' vertex',i3,' is inside body',i3)
	return
6	continue
c
c  call forward routines
	j1=j
	if(grvswt) call pickmg(0,j1,rho(j),ncrn(j),xb1,zb1,
     1 nobs,xg,zg,gc)
	if(magswt)  call pickmg(1,j1,sus(j),ncrn(j),xb1,zb1,
     1 mobs,xh,zh,hc)
100	continue
c
c  gravity errors
	if(.not.grvswt) go to 30
	if(iave.eq.0) go to 19
	ave=0.0
	do 16 i=1,nobs
16	ave=ave+(go(i)-gc(i))
	ave=ave/float(nobs)
	print 17,ave
17	format(' removing',1pg15.3,' mgals from calculated gravity')
	do 18 i=1,nobs
18	gc(i)=gc(i)+ave
19	errt=0.0
	do 20 i=1,nobs
	gerr(i)=go(i)-gc(i)
20	errt=errt+gerr(i)*gerr(i)
	rmsg=sqrt(errt/float(nobs))
c
c  magnetic errors
30	if(.not.magswt) return
	if(jave.eq.0) go to 39
	ave=0.0
	do 36 i=1,mobs
36	ave=ave+(ho(i)-hc(i))
	ave=ave/float(mobs)
	print 37,ave
37	format(' removing',1pg15.3,' nTesla from calculated magnetics')
	do 38 i=1,mobs
38	hc(i)=hc(i)+ave
39	errt=0.0
	do 40 i=1,mobs
	herr(i)=ho(i)-hc(i)
40	errt=errt+herr(i)*herr(i)
	rmsm=sqrt(errt/float(mobs))
	return
	end
c******************************************************************************
	subroutine pickmg(mgswt,ibody,rho,nvert,xb,zb,
     1 nfld,xf,zf,pf)
c  switching routine for the individual forward algorithms
c  rho is the passed linear parameter (density or susceptibility)
c  prism lengths and remanent magnetic values are passed thru common
c  xb,zb is one body coordinate array
c  xf,zf array of locations to be used
	common /model/ dum(500),pl(50,2),rf(50,3),dum1(150)
	common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
	dimension yl(2),vm(3),xb(1),zb(1),xf(1),zf(1),pf(1)
	data d2r/1.745329e-2/
	yl(1)=pl(ibody,1)
	yl(2)=pl(ibody,2)
	if(mgswt.eq.0) go to 10
	rinc=d2r*rf(ibody,2)
	rdec=d2r*(rf(ibody,3)-azmuth)
	rl=cos(rdec)*cos(rinc)
	rm=sin(rdec)*cos(rinc)
	rn=sin(rinc)
c  convert emu/cc to gamma
	rem=rf(ibody,1)*1.e5
c  polarization vector
	vm(1)=rho*ef(1)*ev(1)+rem*rl
	vm(2)=rho*ef(1)*ev(2)+rem*rm
	vm(3)=rho*ef(1)*ev(3)+rem*rn
	ylen=abs(yl(1))
	if(ylen.eq.0.0) ylen=1000.
	call mag2hd(vm,tdir,ylen,nvert,xb,zb,nfld,xf,zf,pf)
	return
10	if(yl(1).eq.0. .and. yl(2).eq.0.) go to 20
	call grv2hd(rho,yl,nvert,xb,zb,nfld,xf,zf,pf)
	return
20	call talwon(rho,nvert,xb,zb,nfld,xf,zf,pf)
	return
	end
c******************************************************************************
	subroutine talwon(rho,nvert,xb,zb,nobs,
     & xf,zf,gc)
c  gravity profile using the Talwani 2-d prism algorithm
c  xb,zb,xf,zf in kilometers, rho in gm/cc
c  reference Talwani et al, JGR, v64, no1.
c  coded by Mike Webring 9/80
	dimension xb(1),zb(1),xf(1),zf(1),gc(1)
	if(rho.eq.0.0) return
	const=2.0*6.673
	pi=3.1415927
	hfpi=pi*.5
	twopi=2.0*pi
	tol=1.e-5
c
	do 200 j=1,nobs
	gtmp=0.0
c
	do 101 i=1,nvert
	capz=0.0
	i1=i+1
	if(i.eq.nvert) i1=1
	x1=xb(i)-xf(j)
	x2=xb(i1)-xf(j)
	z1=zb(i)-zf(j)
	z2=zb(i1)-zf(j)
	ax1=abs(x1)
	ax2=abs(x2)
	if(ax1.lt.tol .and. abs(z1).lt.tol) go to 100
	if(ax2.lt.tol .and. abs(z2).lt.tol) go to 100
	dx=x2-x1
	dz=z2-z1
	t1=atan2(z1,x1)
	t2=atan2(z2,x2)
	if(t1.lt.0.0) t1=t1+twopi
	if(t2.lt.0.0) t2=t2+twopi
	dt=t1-t2
	if(abs(dt).gt.pi) dt=dt-sign(twopi,dt)
	if(abs(dt).lt.tol) go to 100
c
	if(abs(dx).gt.tol) go to 10
	ct2=cos(t2)
	if(abs(ct2).lt.tol) go to 100
	f1=cos(t1)/ct2
	capz=x1*log(f1)
	go to 100
c
10	if(abs(dz).gt.tol) go to 20
	capz=-dt*z1
	go to 100
c
20	tphi=dz/dx
	dxs=dx*dx
	r2=dxs+dz*dz
	if(r2.lt.tol) go to 100
	acs=(dx*dz*x2-dxs*z2)/r2
	if(ax1.gt.tol) go to 30
	f2=sin(t2)-cos(t2)*tphi
	capz=-acs*(-dt+tphi*log(f2))
	go to 100
c
30	if(ax2.gt.tol) go to 40
	f1=sin(t1)-cos(t1)*tphi
	if(f1.le.0.0) then
	print *,j,xf(j),zf(j),xb(i),zb(i),xb(i1),zb(i1)
	stop
	endif
	capz=acs*(dt+tphi*log(f1))
	go to 100
c
40	r1s=x1*x1+z1*z1
	r2s=x2*x2+z2*z2
	capz=acs*(dt+tphi*0.5*log(r2s/r1s))
	go to 100
c
100	gtmp=gtmp+capz
101	continue
200	gc(j)=gc(j)+const*rho*gtmp
	return
	end
c******************************************************************************
	subroutine grv2hd(rho,yin,nvert,xb,zb,nobs,
     & xf,zf,gc)
c  profile gravity calculation of a 2.5d prism
c  yin1 > yin2, z positive down
c  reference Rasmussen and Pedersen, Geophysical Prospecting, v27 p749-760.
c  coded by Mike Webring, 11/81
	dimension x(2),y(2),yin(2),z(2),xzr(2),yr(2,2),f(4,2),ef(8)
	dimension xb(nvert),zb(nvert),xf(nobs),zf(nobs),gc(nobs)
	logical sym
	equivalence (f(1,1),ef(1))
c  rho in gm/cc, distance in kilometer
	const=6.673*rho
	tol=1.e-5
	pi=3.1415926
	twopi=2.*pi
	sym=.false.
	if(yin(1).eq.-yin(2)) sym=.true.
	it=2
	if(sym) it=1
	add=sign(1.,yin(1)*yin(2))
	y(1)=abs(yin(1))
	y(2)=abs(yin(2))
	if(y(1).gt.y(2)) go to 1
	y(1)=abs(yin(2))
	y(2)=abs(yin(1))
1	if(y(1).lt.tol .or. rho.eq.0.0) return
	if(y(2).lt.tol) it=1
c
	do 200 j=1,nobs
	x(1)=xb(1)-xf(j)
	z(1)=zb(1)-zf(j)
	r2=x(1)*x(1)+z(1)*z(1)
	xzr(1)=sqrt(r2)
	yr(1,1)=sqrt(r2+y(1)*y(1))
	yr(1,2)=sqrt(r2+y(2)*y(2))
c
	do 10 k=1,8
10	ef(k)=0.
	do 100 i=1,nvert
	i1=i+1
	if(i.eq.nvert) i1=1
	x(2)=xb(i1)-xf(j)
	z(2)=zb(i1)-zf(j)
	r2=x(2)*x(2)+z(2)*z(2)
	xzr(2)=sqrt(r2)
	dx=x(2)-x(1)
	dz=z(2)-z(1)
	dxz2=dx*dx+dz*dz
	dl=sqrt(dxz2)
	if(dl.lt.tol) go to 90
	cosp=dx/dl
	zn=-cosp
	sinp=dz/dl
	w=cosp*z(1)-sinp*x(1)
	u1=cosp*x(1)+sinp*z(1)
	u2=cosp*x(2)+sinp*z(2)
c
	do 80 ip=1,it
	yr(2,ip)=sqrt(r2+y(ip)*y(ip))
	f(1,ip)=f(1,ip)+zn*y(ip)*alog((u2+yr(2,ip))/(u1+yr(1,ip)))
	f(2,ip)=f(2,ip)+zn*u2*alog((y(ip)+yr(2,ip))/xzr(2))
	f(3,ip)=f(3,ip)+zn*u1*alog((y(ip)+yr(1,ip))/xzr(1))
	xn=u2*y(ip)
	xd=w*yr(2,ip)
	t2=atan2(xn,xd)
	if(t2.lt.0.0) t2=t2+twopi
	xn2=u1*y(ip)
	xd2=w*yr(1,ip)
	t1=atan2(xn2,xd2)
	if(t1.lt.0.0) t1=t1+twopi
	dt=t2-t1
	if(abs(dt).gt.pi) dt=dt-sign(twopi,dt)
	f(4,ip)=f(4,ip)+zn*w*dt
80	continue
90	xzr(1)=xzr(2)
	yr(1,1)=yr(2,1)
	yr(1,2)=yr(2,2)
	x(1)=x(2)
	z(1)=z(2)
100	continue
c
	sum1=f(1,1)-f(4,1)+f(2,1)-f(3,1)
	sum2=sum1
	if(.not.sym) sum2=f(1,2)-f(4,2)+f(2,2)-f(3,2)
	gf=sum1-add*sum2
	gc(j)=gc(j)-const*gf
200	continue
	return
	end
c******************************************************************************
	subroutine mag2hd(vm,tf,y,nvert,xb,zb,nobs,
     & xf,zf,t)
c  computes the total field profile for a 2.5d prism
c  ref: shuey and pasquale, geophysics v38 n3
c  vm is the magnetization vector and tf is measurement unit vector
c  coded by mike webring, 10/81
	complex cdxz,cxz,ct1,ct2,cf,cfln
	dimension vm(3),tf(3),cf(2),x(2),z(2),r(2),
     1 xb(nvert),zb(nvert),xf(nobs),zf(nobs),t(nobs)
	if(y.gt.0.0) go to 1
	print *,' mag2hd: prism length <= zero, setting to 1000'
	y=1000.
1	ys=y*y
	do 50 iob=1,nobs
	px=0.
	pz=0.
	q=0.
	x(1)=xb(1)-xf(iob)
	z(1)=zb(1)-zf(iob)
	r(1)=sqrt(x(1)*x(1)+ys+z(1)*z(1))
c
	do 20 j=1,nvert
	j1=j+1
	if(j.eq.nvert) j1=1
	x(2)=xb(j1)-xf(iob)
	z(2)=zb(j1)-zf(iob)
	  if(x(2).eq.0.0 .and. z(2).eq.0.0) then
	  print *,' mag2hd: observation coincident with prism corner'
	  return
	endif
	r(2)=sqrt(x(2)*x(2)+ys+z(2)*z(2))
	dx=x(2)-x(1)
	dz=z(2)-z(1)
	if(dx.eq.0.0 .and. dz.eq.0.0) then
	  print *,' mag2hd: zero length segment, vertex ',j
	  return
	endif
	cdxz=cmplx(dx,dz)
	 do 10 k=1,2
	 cxz=cmplx(x(k),z(k))
	 t1=(1.0+r(k)/y)
	 ct1=cmplx(t1,0.)
	 t2=(x(k)*dz-z(k)*dx)/ys
	 ct2=cmplx(0.,t2)
10	 cf(k)=(cdxz*ct1)/cxz+ct2
	cfln=(clog(cf(2)/cf(1)))/cdxz
	px=px+dz*real(cfln)
	pz=pz-dx*aimag(cfln)
	q=q-dx*real(cfln)
	x(1)=x(2)
	z(1)=z(2)
20	r(1)=r(2)
c
	hx=2.0*(vm(1)*px+vm(3)*q)
	hy=2.0*vm(2)*(pz-px)
	hz=2.0*(vm(1)*q-vm(3)*pz)
50	t(iob)=t(iob)+hx*tf(1)+hy*tf(2)+hz*tf(3)
	return
	end
c******************************************************************************
	subroutine nlctl(ng,garr,nh,harr)
c  setup routine for simple parameters: x,z,rho,sus
c  watch place holding ng,nh vs. real ngrv,nmag dimensions
	common /parm1/  nvert,nbody,nobs1,mobs1,grvswt,magswt,
     1               iave(2),vex,istatn,iparsh,labelv,numbod
	common /model/  xb(200),zb(200),rho(50),sus(50),dum(400)
	common /swork/  a(8000),err(400),s(20)
	common /freel/  weight,rms1,iter,xk1,nfreex,nfreez,nfreer,nfrees,
     1               ifree(20),jfree(20),kfree(20),lfree(20)
	dimension garr(ng,5),harr(nh,5)
	logical grvswt,magswt
	data delta/1.0e-3/
c
c  get free parameter lists
	call getpar(npart,ier)
	if(ier.eq.1) return
c
c  get limits for active profile segment
100	call getact(igs,ngrv,ihs,nmag)
	no=ngrv+nmag
	if(no.le.2) return
c
c  get error vector
	do 101 i=1,ng
101	garr(i,4)=0.0
	do 102 i=1,nh
102	harr(i,4)=0.0
	iad=1
	if(grvswt)  then
	            call mardvr(0,xb,zb,rho,igs,ngrv,
     1           garr(igs,4),err(iad),rmsg)
	            iad=ngrv+1
	endif
	if(magswt)  call mardvr(1,xb,zb,sus,ihs,nmag,
     1           harr(ihs,4),err(iad),rmsh)
c
c  setup jacobian
	iad=1
	do 105 i=1,no*npart
105	a(i)=0.0
c
	if( nfreex.lt.1 ) go to 115
	dx=delta
	dz=0.0
	do 110 i=1,nfreex
	ivert=ifree(i)
	if(grvswt) call parshl(0,dx,dz,ivert,
     1          ngrv,garr(igs,1),garr(igs,2),a(iad))
	iad=iad+ngrv
	if(magswt) call parshl(1,dx,dz,ivert,
     1          nmag,harr(ihs,1),harr(ihs,2),a(iad))
110	iad=iad+nmag
c
115	if( nfreez.lt.1 ) go to 130
	dx=0.0
	dz=delta
	do 120 i=1,nfreez
	ivert=jfree(i)
	if(grvswt) call parshl(0,dx,dz,ivert,
     1          ngrv,garr(igs,1),garr(igs,2),a(iad))
	iad=iad+ngrv
	if(magswt) call parshl(1,dx,dz,ivert,
     1          nmag,harr(ihs,1),harr(ihs,2),a(iad))
120	iad=iad+nmag
c
130	dx=0.0
	dz=0.0
	if(nfreer.eq.0) go to 150
	do 140 i=1,nfreer
	iparm=kfree(i)
	if(grvswt) call parshl(0,dx,dz,iparm,
     1          ngrv,garr(igs,1),garr(igs,2),a(iad))
	iad=iad+ngrv
	if(magswt) then
	do 135 k=1,nmag
	a(iad)=0.0
135	iad=iad+1
	endif
140	continue
c
150	if(nfrees.eq.0) go to 170
	do 160 i=1,nfrees
	if(grvswt) then
	do 155 k=1,ngrv
	a(iad)=0.0
155	iad=iad+1
	endif
	iparm=lfree(i)
	call parshl(1,dx,dz,iparm,
     1    nmag,harr(ihs,1),harr(ihs,2),a(iad))
	iad=iad+nmag
160	continue
c
170	if(iparsh.eq.1) then
	  print 171
171	  format(' plot partial derivatives ? '$)
	  if(noyes(l).eq.1) call advr(igs,ngrv,ihs,nmag,no,npart,a,err)
	endif
c
c  relative weighting of gravity to magnetics
	npxz=nfreex+nfreez
	call relsca(weight,ngrv,nmag,npxz,no,npart,a,err)
	if(iparsh.eq.1) then
	  print 172
172	  format(' plot weighted partial derivatives ? '$)
	  if(noyes(l).eq.1) call advr(igs,ngrv,ihs,nmag,no,npart,a,err)
	endif
c
	if(iter.eq.1) then
	if(grvswt .and. magswt) then
	  rms1=(weight*rmsg+rmsh)/(weight+1.0)
	  print 175,rms1
175	  format(1x,' weighted starting rms error',1pg15.5)
	else
	  if(grvswt) rms1=rmsg
	  if(magswt) rms1=rmsh
	  print 176,rms1
176	  format(' starting rms error',1pg15.5)
	endif
	endif
c
c  scale jacobian
	call cscale(no,npart,a,s)
	if(iparsh.eq.1) then
	  print 173
173	  format(' plot scaled system ? '$)
	  if(noyes(l).eq.1) call advr(igs,ngrv,ihs,nmag,no,npart,a,err)
	endif
c
	if(iter.gt.1) go to 200
	print 180
180	format(' print parameter correlations ? '$)
	if(noyes(l).eq.1) then
	  call crlmtx(6,no,npart,a)
	  print 190
190	  format(' continue with the inversion ? '$)
	  if(noyes(l).eq.0) return
	endif
c
c  generate new parameters
200	call marq(no,npart,igs,ngrv,ihs,nmag)
	iter=iter+1
	print 210
210	format(' another iteration ? '$)
	if(noyes(l).eq.1) go to 100
999	return
	end
c******************************************************************************
	subroutine getact(igs,ngrv,ihs,nmag)
        common /field/  nf(10),f(5120)
	common /model/  xb(200),zb(200),rho(50),sus(50),dum(400)
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,dum1(7)
	common /vertx/  icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	common /freel/  weight,rms1,iter,xk1,nfreex,nfreez,nfreer,nfrees,
     1               ifree(20),jfree(20),kfree(20),lfree(20)
	logical grvswt,magswt
c  the partial derivatives should be approaching zero at pfact*depth.
	pfact=5.
	xmax=-1.e5
	xmin= 1.e5
	avedat= 1.e5
	if(grvswt) then
	  aveg=0.0
	  do 5 i=0,nobs-1
5	  aveg=aveg+f(nf(2)+i)
	  avedat=aveg/nobs
	endif
	if(magswt) then
	  aveh=0.0
	  do 10 i=0,mobs-1
10	  aveh=aveh+f(nf(7)+i)
	  aveh=aveh/mobs
	  avedat=amin1(avedat,aveh)
	endif
c  free x's
	do 20 ix=1,nfreex
	x=xb(ifree(ix))
	depth=abs(zb(ifree(ix))-avedat)
	x1=x-pfact*depth
	x2=x+pfact*depth
	if(x1.lt.xmin) xmin=x1
	if(x2.gt.xmax) xmax=x2
20      continue 
c  free z's
	do 30 iz=1,nfreez
	x=xb(jfree(iz))
	depth=abs(zb(jfree(iz))-avedat)
	x1=x-pfact*depth
	x2=x+pfact*depth
	if(x1.lt.xmin) xmin=x1
	if(x2.gt.xmax) xmax=x2
30      continue 
c  free rho and susceptibility
	do 60 k=1,2
	 nfree=0
	 if(k.eq.1 .and. grvswt) nfree=nfreer
	 if(k.eq.2 .and. magswt) nfree=nfrees
	do 50 j=1,nfree
	 ibody=kfree(j)	
	 if(k.eq.2) ibody=lfree(j)
	do 40 i=1,ncrn(ibody)
	iadr=icrn(loca(ibody)+i)
	x=xb(iadr)
	depth=abs(zb(iadr)-avedat)
	x1=x-pfact*depth
	x2=x+pfact*depth
	if(x1.lt.xmin) xmin=x1
	if(x2.gt.xmax) xmax=x2
40      continue
50      continue
60      continue 
c
c  appropriate indices
	igs=1
	ihs=1
	ngrv=0
	nmag=0
	if(grvswt) then
	 ige=nobs
	 x1=f(nf(1))
	do 70 i=1,nobs-1
	 x2=f(nf(1)+i)
	 if(xmin.ge.x1 .and. xmin.lt.x2) igs=i
	 if(xmax.gt.x1 .and. xmax.le.x2) ige=i+1
	 x1=x2
70       continue 
	 ngrv=ige-igs+1
	 xs = f(nf(1)+igs)
	 xe = f(nf(1)+ige-1)
	 print *,' active  gravity range (km)',xs,xe
	endif
c
	if(magswt) then
	 ihe=mobs
	 x1=f(nf(6))
	do 80 i=1,mobs-1
	 x2=f(nf(6)+i)
	 if(xmin.ge.x1 .and. xmin.lt.x2) ihs=i
	 if(xmax.gt.x1 .and. xmax.le.x2) ihe=i+1
	 x1=x2
80       continue 
	 nmag=ihe-ihs+1
	 xs = f(nf(6)+ihs)
	 xe = f(nf(6)+ihe-1)
	 print *,' active magnetic range (km)',xs,xe
	endif
	return
	end
c******************************************************************************
	subroutine getpar(npart,ier)
c  enter the free parameters x,z,rho,sus
	common /parm1/  nvert,nbody,nobs,mobs,grvswt,magswt,dum(7)
	common /freel/  weight,rms1,iter,xk1,nfreex,nfreez,nfreer,nfrees,
     1               ifree(20),jfree(20),kfree(20),lfree(20)
	logical grvswt,magswt
	data itty/5/
c
	ier=0
5	npart=nfreex+nfreez+nfreer+nfrees
	if(npart.ne.0) then
	  print 6
6	  format(' reset free parameters ? '$)
	  if(noyes(l).eq.0) go to 200
	endif
	npart=0
	nfreer=0
	nfrees=0
	iter=1
9	print  10
10	format(' enter number of free x vertices:'$)
	read(itty,*,end=999,err=999) nfreex
	if(nfreex.lt.0) go to 101
	ntmp=nfreex
	if(nfreex.gt.nvert .or. ntmp.gt.20) go to 120
	if(nfreex.eq.0) go to 15
	print  24
	read(itty,*,end=9,err=9) (ifree(i),i=1,nfreex)
	do 12 i=1,nfreex
	if(ifree(i).lt.1 .or. ifree(i).gt.nvert) go to 9
12	continue
c
15	print  20
20	format(' enter number of free z vertices:'$)
	read(itty,*,end=999,err=999) nfreez
	if(nfreez.lt.0) go to 101
	ntmp=nfreez+nfreex
	if(nfreez.gt.nvert .or. ntmp.gt.20) go to 120
	if(nfreez.eq.0) go to 25
	print  24
24	format(' enter indices : '$)
	read(itty,*,end=15,err=15) (jfree(i),i=1,nfreez)
	do 22 i=1,nfreez
	if(jfree(i).lt.1 .or. jfree(i).gt.nvert) go to 15
22	continue
c
25	if(.not.grvswt) go to 70
	print  50
50	format(' enter number of free densities :'$)
	read(itty,*,end=999,err=999) nfreer
	if(nfreer.lt.0) go to 101
	ntmp=nfreer+nfreez+nfreex
	if(nfreer.gt.nbody .or. ntmp.gt.20) go to 120
	if(nfreer.eq.0) go to 70
	print  24
	read(itty,*,end=25,err=25) (kfree(i),i=1,nfreer)
	do 60 i=1,nfreer
	if(kfree(i).lt.1 .or. kfree(i).gt.nbody) go to 25
60	continue
c
70	if(.not.magswt) go to 100
	print  80
80	format(' enter number of free susceptibilites :'$)
	read(itty,*,end=999,err=999) nfrees
	if(nfrees.lt.0) go to 101
	ntmp=nfrees+nfreer+nfreez+nfreex
	if(nfrees.gt.nbody .or. ntmp.gt.20) go to 120
	if(nfrees.eq.0) go to 100
	print  24
	read(itty,*,end=70,err=70) (lfree(i),i=1,nfrees)
	do 90 i=1,nfrees
	if(lfree(i).lt.1 .or. lfree(i).gt.nbody) go to 70
90	continue
c
c  npartials check
100	npart=nfreex+nfreez+nfreer+nfrees
101	if(npart.gt.0) go to 110
	print *,' none set'
	nfreex=0
	nfreez=0
	nfreer=0
	nfrees=0
	go to 999
110	if(npart.le.20) go to 200
	print *,' maximum number free is 20'
	go to 999
120	print *,' exceeded either number of vertices or bodies' 
	print *,' or a total of 20 parameters free.'
	go to 999
c
200	return
999	ier=1
	return
	end
c******************************************************************************
	subroutine parshl(mgswt,dx,dz,iparm,nfld,xf,zf,pf)
c  partial derivatives for x,z,rho,sus via central differences
	common /parm1/  nvert,nbody,dum1(11)
	common /model/  xb(200),zb(200),rho(50),sus(50),pl(50,2),
     1 rem(50,3),dum(150)
	common /vertx/  icrn(500),ncrn(150),loca(150),dum2(201)
	common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
        dimension xb1(200),zb1(200),tmp(400)
	dimension xf(nfld),zf(nfld),pf(nfld)
c        data d2r/1.745329e-2/
	do 1 i=1,nfld
	tmp(i)=0.0
1	pf(i)=0.0
c
	if(dx.eq.0.0 .and. dz.eq.0.0) go to 200
	do 100 jbody=1,nbody
	is=1
	ie=ncrn(jbody)
	if(iparm.le.0) go to 21
	do 10 ip=1,ncrn(jbody)
	if(icrn(loca(jbody)+ip).eq.iparm) go to 20
10	continue
	go to 100
20	is=ip
	ie=ip
21	if(mgswt.eq.0) then
	  alin=rho(jbody)
	else
	  alin=sus(jbody)
	endif
	do 30 i=1,ncrn(jbody)
	xb1(i)=xb(icrn(loca(jbody)+i))
30	zb1(i)=zb(icrn(loca(jbody)+i))
	do 40 i=is,ie
	xb1(i)=xb1(i)+dx
40	zb1(i)=zb1(i)+dz
	j1=jbody
	call pickmg(mgswt,j1,alin,ncrn(jbody),xb1,zb1,
     1 nfld,xf,zf,pf)
	do 50 i=is,ie
	xb1(i)=xb1(i)-2.0*dx
50	zb1(i)=zb1(i)-2.0*dz
	call pickmg(mgswt,j1,alin,ncrn(jbody),xb1,zb1,
     1 nfld,xf,zf,tmp)
100	continue
	delta=dx
	if(dx.eq.0.0) delta=dz
	delta=1.0/(2.0*abs(delta))
	go to 300
c
200	do 250 i=1,ncrn(iparm)
	  xb1(i)=xb(icrn(loca(iparm)+i))
	  zb1(i)=zb(icrn(loca(iparm)+i))
250     continue 
	if(mgswt.eq.0) then
	  alin=1.0
	  call pickmg(0,iparm,alin,ncrn(iparm),xb1,zb1,nfld,xf,zf,pf)
	  return
	else
	  h=abs(sus(iparm)*ef(1)) + abs(rem(iparm,1)*1.0e5)
	  ds=h*1.0e-5
	  if(ds.lt.1.e-10) ds=1.e-10
	  alin=sus(iparm)+ds
	  call pickmg(1,iparm,alin,ncrn(iparm),xb1,zb1,nfld,xf,zf,pf)
	  alin=sus(iparm)-ds
	  call pickmg(1,iparm,alin,ncrn(iparm),xb1,zb1,nfld,xf,zf,tmp)
	  delta=1.0/(2.0*ds)
	endif
c
300	do 350 i=1,nfld
	pf(i)=delta*(pf(i)-tmp(i))
350     continue
	return
	end
c******************************************************************************
	subroutine crlmtx(jdev,m,n,a)
c  print the triangular correlation matrix
	dimension icl(20),a(m,n)
	write(jdev,5)
5	format(1x,' a=')
	do 20 j=2,n
	do 10 i=1,j-1
10	icl(i)=int( 100.0*correl(m,a(1,j),a(1,i)) + 0.5)
	write(jdev,15) j,(icl(i),i=1,j-1)
15	format(1x,i3,20i3)
20	continue
	write(jdev,30)
30	format(/,1x,' b=  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15')
	return
	end
c******************************************************************************
	subroutine relsca(weight,ngrv,nmag,npxz,no,np,a,err)
c  row scaling for relative mag/grv weight
	common /parm1/ fill(4),grvswt,magswt,dum(7)
	dimension a(no,np),err(no)
	logical grvswt,magswt
	weight=1.0
	if(.not.(grvswt.and.magswt)) return
	call ammi(ngrv,err,        gmin,gmax,mn,mx,1)
	call ammi(nmag,err(ngrv+1),hmin,hmax,mn,mx,1)
	gmax=amax1(abs(gmin),abs(gmax))
	hmax=amax1(abs(hmin),abs(hmax))
	ratio=1.0
	if(gmax*hmax.gt.1.e-10) ratio=gmax/hmax
	do 10 j=1,npxz
	  call ammi(ngrv,a(1,j),     gmin,gmax,mn,mx,1)
	  call ammi(nmag,a(ngrv+1,j),hmin,hmax,mn,mx,1)
	  gmax=amax1(abs(gmin),abs(gmax))
	  hmax=amax1(abs(hmin),abs(hmax))
	  r=1.0
	  if(gmax*hmax.gt.1.e-10) r=gmax/hmax
	  ratio=ratio+r
10      continue 
	ratio=ratio/float(npxz+1)
c
	relwt=weight*ratio
	do 50 j=1,np
	do 50 i=ngrv+1,no
50	a(i,j)=relwt*a(i,j)
	do 100 i=ngrv+1,no
100	err(i)=relwt*err(i)
	return
	end
c******************************************************************************
	subroutine cscale(m,n,a,s)
c  column scaling for matrix a
	dimension s(n),a(1)
	tol=1.e-10
	do 10 j=1,n
	k=(j-1)*m+1
	pmin=a(k)
	pmax=a(k)
	do 5 i=1,m
	if(a(k).gt.pmax) pmax=a(k)
	if(a(k).lt.pmin) pmin=a(k)
5	k=k+1
	s(j)=amax1(abs(pmax),abs(pmin))
	if(s(j).lt.tol) s(j)=1.0
10	s(j)=1.0/s(j)
	k=1
	do 20 j=1,n
	do 20 i=1,m
	a(k)=a(k)*s(j)
20	k=k+1
	return
	end
c******************************************************************************
	subroutine marq(no,np,igs,ngrv,ihs,nmag)
        common /field/  nf(10),f(5120)
	common /parm1/  nv,nb,ifill(2),grvswt,magswt,dum(7)
	common /model/  xb(200),zb(200),rho(50),sus(50),dum1(400)
	common /freel/  weight,rms1,iter,xk1,nfx,nfz,nfr,nfs,
     1               ifx(20),ifz(20),ifr(20),ifs(20)
	common /swork/  a(8000),err(400),s(20)
       common /swork2/ u(8000),v(400)
c
	dimension q(20),  work(80),ute(20),qute(20),dp(20),dg(20)
	equivalence (work,ute),    (work(21),qute),
     1           (work(41),dp), (work(61),dg)
	dimension xb1(200),zb1(200),rho1(50),sus1(50)
	dimension nf1(10),xk(3),rms(3),ang(3),dp2(3)
c
	character cp*1
	logical grvswt,magswt
	data attn/4./,mxridg/10/
c  eigenvalue attenuation rate set by attn
c
	do 5 i=1,5
	  nf1(i)  = nf(i)  +(igs-1)
	  nf1(i+5)= nf(i+5)+(ihs-1)
5       continue
	icount=0
	istop=0
	limit=0
	init=1
	ii=3
c
	call svd(a,no,np,q,u,v,work)
	if(q(1).eq.0.0) then
	  print *,' no independent parameters'
	  return
	endif
	print *,' scaled eigenvalues'
	do 7 i=1,np
	  ute(i)=q(i)/q(1)
7       continue
	print 10,(ute(i),i=1,np)
10	format(4(1x,5f8.4,/))
c
	call atb(u,no,np,err,no,1,ute)
	call ata(u,no,np,a)
	call chki(a,np,icode)
	call ata(v,np,np,a)
	call chki(a,np,jcode)
	if(icode.lt.6) print *,'svd: UtU',icode
	if(jcode.lt.6) print *,'svd: VtV',jcode
c
	call ammi(np,q,qmin,qmax,mn,mx,1)
	xk(3) = 1.5*qmax
	if (iter.gt.1)  xk(3) = 1.5*xk1
c
c  steepest descent vector
	damp=10.*qmax**attn
	do 50 i=1,np
50	qute(i)=ute(i)/(q(i)+damp)
	call axb(v,np,np,qute,np,1,dg)
	dgl=0.
	do 80 i=1,np
	dg(i)=dg(i)*s(i)
80	dgl=dgl+dg(i)*dg(i)
	dgl=sqrt(dgl)
c
c  system damping loop
	print  90
90	format('  cutoff     rms      angle     delta p')
100	icount=icount+1
	if((limit.eq.1  .and. init.eq.0) .or. icount.gt.mxridg) then
	  istop=1
	  call ammi(3,rms,emin,emax,ii,mx,1)
	  call ammi(3,ang,amin,amax,iia,mx,1)
	  if(iia.gt.ii) ii=iia
	  if (rms(ii).gt.rms1) ii=3
	  thres=xk(ii)/q(1)
	  print 101,thres
101	  format(' choosing eigenvalue cutoff',f8.4)
	endif
c
c  get new correction vector
	damp=xk(ii)**attn
	do 104 i=1,np
	q2=q(i)**(attn-1.)
104	qute(i)=ute(i)*(q2/(q2*q(i)+damp))
	call axb(v,np,np,qute,np,1,dp)
	do 105 i=1,np
105	dp(i)=dp(i)*s(i)
c
c  assemble and check parameter vector
        call corvec(limit,dp,xb1,zb1,rho1,sus1,istop)
c
c  forward calculation
200	if(grvswt) call mardvr(0,xb1,zb1,rho1,igs,ngrv,
     1          f(nf1(4)),f(nf1(5)),rmsg)
	if(magswt) call mardvr(1,xb1,zb1,sus1,ihs,nmag,
     1          f(nf1(9)),f(nf1(10)),rmsm)
	rms(ii)=rmsg
	if(magswt) rms(ii)=rmsm
	if(grvswt .and. magswt) rms(ii)=(weight*rmsg+rmsm)/(weight+1.0)
	dp2(ii)=0.0
	an=0.0
	dpl=0.0
	do 210 i=1,np
	  if(i.le.nfx+nfz) dp2(ii) = dp2(ii) + dp(i)*dp(i)
	  dpl = dpl + dp(i)*dp(i)
210	  an  = an  + dp(i)*dg(i)
	dpl = sqrt(dpl)
	if(dpl.lt.1.e-10) go to 300
	ad = dpl*dgl
	ang(ii) = acosd2(an,ad)
	if(nfx+nfz.gt.0) dp2(ii) = sqrt( dp2(ii)/float(nfx+nfz) )
	thres = xk(ii)/q(1)
	print  220, thres, rms(ii), ang(ii), dp2(ii)
220	format(1x,f8.4,g12.4,f7.2,f10.3)
	if(istop.eq.1) go to 300
c
c  iteration control
	if(init.eq.1 .and. ii.gt.1) then
	  ii=ii-1
	  xk(ii)=xk(ii+1)/2.
	else
	  init=0
	  call bictl(xk,rms,ang,dp2,ii,ireset)
	  if(ireset.eq.-1) limit=1
	endif
	go to 100
c
c  print new parameter values
300	print  310
310	format(' parm   change       new value')
	do 390 i=1,np
	if(i.gt.nfx) go to 320
	cp='x'
	j=ifx(i)
	p=xb1(j)
	go to 380
320	if(i.gt.(nfx+nfz)) go to 330
	cp='z'
	j=ifz(i-nfx)
	p=zb1(j)
	go to 380
330	if(i.gt.(nfx+nfz+nfr)) go to 340
	cp='r'
	j=ifr(i-(nfx+nfz))
	p=rho1(j)
	go to 380
340	cp='s'
	j=ifs(i-(nfx+nfz+nfr))
	p=sus1(j)
380	print  381,cp,j,dp(i),p
381	format(1x,a1,i3,5x,1p2(g12.4,1x))
390	continue
c
c  save new parameter values
	rmsimp=0.0
	if(rms1.gt.0.0) rmsimp=100.*(rms1-rms(ii))/rms1
	print  410,iter,rms(ii),rmsimp
410	format(' iteration ',i2,',  rms error is  ',1pg15.4,/,
     1 ' percent improvement is  ',1pg15.4)
	print  420
420	format(' save new parameter values ?'$)
	if(noyes(l).eq.0) return
	do 430 i=1,nv
	xb(i)=xb1(i)
430	zb(i)=zb1(i)
	do 440 i=1,nb
	if(magswt) sus(i)=sus1(i)
	if(grvswt) rho(i)=rho1(i)
440	continue
	xk1=xk(ii)
	rms1=rms(ii)
	return
	end
c******************************************************************************
        subroutine corvec(limit,dp,xb1,zb1,rho1,sus1,iprnt)
        common /field/  nf(10),f(5120)
	common /parm1/  nv,nb,ifill(2),grvswt,magswt,dum(7)
	common /model/  xb(200),zb(200),rho(50),sus(50),dum1(400)
	common /freel/  fill(4),nfx,nfz,nfr,nfs,
     1               ifx(20),ifz(20),ifr(20),ifs(20)
c
	dimension dp(20),xb1(200),zb1(200),rho1(50),sus1(50)
	logical grvswt,magswt,limitp
	data rlimit/.05/, xfrac,zfrac/50.0,10.0/
c
	limit=0
c  get x&z excursion limits
	x1=0.0
	x2=0.0
	if(grvswt) x1 = abs(f(nf(2)-1)-f(nf(1)) )
	if(magswt) x2 = abs(f(nf(7)-1)-f(nf(6)) )	
	xlimit = amax1(x1,x2) / xfrac
	call ammi(nv,zb,zmn,zmx,mn,mx,1)
	zlimit = (zmx-zmn) / zfrac
c
c  assemble parameter vector
	do 10 i=1,nv
	xb1(i)=xb(i)
10	zb1(i)=zb(i)
	do 20 i=1,nb
	sus1(i)=sus(i)
20	rho1(i)=rho(i)
c
c  add in correction vector and check excursion
	ip=0
c  x parameter
	do 75 i=1,nfx
	limitx=0
	ip=ip+1
	j=ifx(ip)
50	  absdp = abs(dp(ip))
	  xb1(j)= xb(j)+dp(ip)
	  if(absdp.gt.xlimit .or. limitp(j,nb,xb1,zb1,0)) then
	    dp(ip)=dp(ip)/2.0
	    if(absdp.lt.0.01) dp(ip)=0.0
	    limitx=j
	    dpx=dp(ip)
	  if(dpx.ne.0.0) go to 50
	endif
	if(limitx.gt.0) then
	  if(iprnt.eq.1) 
     1   print *,' limited x excursion of parameter',limitx,' to',dpx
	  limit=1
	  endif
75      continue 
c  z parameter
	do 150 i=1,nfz
	limitz=0
	ip=ip+1
	j=ifz(i)
100	  absdp=abs(dp(ip))
	  zb1(j)=zb(j)+dp(ip)
	  if(absdp.gt.zlimit .or. limitp(j,nb,xb1,zb1,0)) then
	    dp(ip)=dp(ip)/2.0
	    if(absdp.lt.0.01) dp(ip)=0.0
	    limitz=j
	    dpz=dp(ip)
	if(dp(ip).ne.0.0) go to 100
	endif
	if(limitz.gt.0) then
	  if(iprnt.eq.1) 
     1  print *,' limited z excursion of parameter',limitz,' to',dpz
	  limit=1
	  endif
150     continue 
c  rho
	do 200 i=1,nfr
	  ip=ip+1
	  j=ifr(i)
	  if(abs(dp(ip)).gt.rlimit) dp(ip)=sign(rlimit,dp(ip))
	  rho1(j)=rho(j)+dp(ip)
200     continue
c  susceptibility
	do 300 i=1,nfs
	  ip=ip+1
	  j=ifs(i)
	  sus1(j)=sus(j)+dp(ip)
300     continue 
c
	return
	end
c******************************************************************************
	subroutine mardvr(mgswt,xb,zb,alin,istart,np,pf,pe,rms)
c  forward driver for the inversion section where xb,zb,alin are 
c  updated parameter vectors.
        common /field/ nf(10),f(5120)
	common /parm1/ nvert,nbody,nobs,mobs,grvswt,magswt,iave(2),
     1 dum(5)
	common /vertx/ icrn(500),ncrn(150),loca(150),nbord,ibord(200)
	dimension xb(200),zb(200),xb1(200),zb1(200),alin(50)
	dimension pf(np),pe(np)
        logical grvswt,magswt
c  initialize output
	rms=1.e30
	do 10 i=1,np
10	pf(i)=0.0
c  find observed data
	noff=0
	if(mgswt.eq.1) noff=5
	ix = nf(1+noff) + (istart-1)
	iz = nf(2+noff) + (istart-1)
	iobs=nf(3+noff) + (istart-1)
c
c  call forward selection routine
	do 100 j=1,nbody
	do 50 i=1,ncrn(j)
	iadr=icrn(loca(j)+i)
	xb1(i)=xb(iadr)
50	zb1(i)=zb(iadr)
	jbody=j
	call pickmg(mgswt,jbody,alin(j),ncrn(j),xb1,zb1,
     1  np,f(ix),f(iz),pf)
100	continue
c
c  float datum
	if(iave(mgswt+1).eq.0) go to 150
	ave=0.0
	i1=iobs
	do 110 i=1,np
	ave=ave+f(i1)-pf(i)
110	i1=i1+1
	ave=ave/float(np)
	do 120 i=1,np
120	pf(i)=pf(i)+ave
c
c  error vector
150	errt=0.0
	i1=iobs
	do 160 i=1,np
	pe(i)=f(i1)-pf(i)
	errt=errt+pe(i)*pe(i)
160	i1=i1+1
	rms=sqrt(errt/float(np))
	return
	end
c******************************************************************************
	subroutine bictl(x,err,ang,dp2,ii,ireset)
c  locate the next eigenvalue cutoff by checking the rms error and the 
c  angle formed between test vectors and the steepest descent vector 
c  and choose the more conservative minimum.
c  ireset=-1 stop,             ireset=0 interval division,
c  ireset=1 interval addition, ireset=2 balance intervals
	dimension x(3),err(3),ang(3),dp2(3)
	ireset=-1
	if(x(3).le.x(2) .or. x(2).le.x(1)) return
c  check for convergence
	if(err(2).eq.0.0) return
	derr = ( abs(err(3)-err(2)) + abs(err(2)-err(1)) )/ err(2)
	if(derr.lt.0.001) return
c  balance intervals
	ratio=(x(3)-x(2))/(x(2)-x(1))
	if(ratio.gt.5.0 .or. ratio.lt.0.2) then
	  ireset=2
	  ii=2
	  x(2)=0.5*(x(3)+x(1))
	  go to 100
	endif
c
	icase=0
	se2=(err(2)-err(1))/(x(2)-x(1))
	se3=(err(3)-err(2))/(x(3)-x(2))
	if(se2.gt.0. .and. se3.gt.0.0) icase=1
	if(se2.lt.0. .and. se3.lt.0.0) icase=4
	if(icase.eq.0) then
	  icase=2
	  if(abs(se2).gt.se3) icase=3
	endif
c
	jcase=0
	sa2=(ang(2)-ang(1))/(x(2)-x(1))
	sa3=(ang(3)-ang(2))/(x(3)-x(2))
	if(sa2.gt.0. .and. sa3.gt.0.0) jcase=1
	if(sa2.lt.0. .and. sa3.lt.0.0) jcase=4
	if(jcase.eq.0) then
	  jcase=2
	  if(abs(sa2).gt.sa3) jcase=3
	endif
c
c  override an interval division by a good angle improvement
	if (icase.ne.4) then
	  dangl = abs(ang(3)-ang(2)) + abs(ang(2)-ang(1))
	  if (icase.le.2 .and. (jcase.gt.2 .and. dangl.gt.5.0) )
     1  icase=(icase+jcase)/2
	endif
c
	ireset=1
	if(icase.eq.1) then
	  x(3)=x(2)
	  x(2)=x(1)
	  x(1)=x(2)/2.
	  err(3)=err(2)
	  err(2)=err(1)
	  ang(3)=ang(2)
	  ang(2)=ang(1)
	  dp2(3)=dp2(2)
	  dp2(2)=dp2(1)
	else if(icase.eq.2 .or. icase.eq.3) then
	  if(icase.eq.2) x(1)=(x(1)+x(2))/2.
	  if(icase.eq.3) x(3)=(x(2)+x(3))/2.
	  ireset=0
	else if(icase.eq.4) then
	  x(1)=x(2)
	  x(2)=x(3)
	  x(3)=1.5*(x(2)-x(1))+x(2)
	  err(1)=err(2)
	  err(2)=err(3)
	  ang(1)=ang(2)
	  ang(2)=ang(3)
	  dp2(1)=dp2(2)
	  dp2(2)=dp2(3)
	endif
	ii=1
	if(icase.gt.2) ii=3
100	err(ii)=1.0e30
	ang(ii)=90.0
	return
	end
c******************************************************************************
	logical function limitp(m,nbody,xb,zb,iprnt)
	common /vertx/ lookup(500),ncrn(150),loca(150),dum(201)
	dimension lint(2),xb(1),zb(1)
	logical planar
	limitp=.false.
	do 10 j=1,nbody
	  do 10 i=1,ncrn(j)
	  k=lookup(loca(j)+i)
	  if(k.eq.m) then
	    ic=loca(j)+1
	    limitp=.not.planar(ncrn(j),xb,zb,lookup(ic),lint)
	    if(iprnt.ne.0) then
	      print 20, j,lint
	    endif
	    if(limitp) return
	  endif
10	continue
20	format(' body',i3,' has intersecting line segments',i3,'and',i3)
	return
	end
c******************************************************************************
	logical function planar(n,x,z,lookup,lint)
c  check polygon for intersecting line segments
	dimension lint(2),lookup(n),x(1),z(1)
	lint(1)=0
	lint(2)=0
	planar=.true.
	if(n.eq.3) return
c        tolfac=.0001
        tolfac=.000001
c  first line segment
	do 10 line1=1,n-2
	i2=line1+1
	if(i2.gt.n) i2=1
	i=lookup(line1)
	i2=lookup(i2)
	dz=z(i2)-z(i)
	dx=x(i2)-x(i)
	adx=abs(dx)
	tol=tolfac*(adx+abs(dz)) + 1.e-10
	if(adx.lt.tol) dx=tol
	am=dz/dx
	amn=amin1(x(i),x(i2))-tol
	amx=amax1(x(i),x(i2))+tol
	aymn=amin1(z(i),z(i2))-tol
	aymx=amax1(z(i),z(i2))+tol
	k=n
	if(line1.eq.1) k=n-1
c  compared to segments: +2 to n segment-1
	do 10 line2=line1+2,k
	j2=line2+1
	if(j2.gt.n) j2=1
	j=lookup(line2)
	j2=lookup(j2)
	dz=z(j2)-z(j)
	dx=x(j2)-x(j)
	adx=abs(dx)
	tol=tolfac*(adx+abs(dz)) + 1.e-10
	if(adx.lt.tol) dx=tol
	bm=dz/dx
	dslope=am-bm
	if(abs(dslope).lt.tolfac) go to 10
	rs=1.0/dslope
	xi=(rs*am)*x(i)-(rs*bm)*x(j)+rs*(z(j)-z(i))
	yi=am*xi-am*x(i)+z(i)
	bmn=amin1(x(j),x(j2))-tol
	bmx=amax1(x(j),x(j2))+tol
	bymn=amin1(z(j),z(j2))-tol
	bymx=amax1(z(j),z(j2))+tol
c  test for intersection within segments
	if(xi.lt.amn .or. xi.lt.bmn) go to 10
	if(xi.gt.amx .or. xi.gt.bmx) go to 10
	if(yi.lt.aymn .or. yi.lt.bymn) go to 10
	if(yi.gt.aymx .or. yi.gt.bymx) go to 10
c  intersection between line segments
	lint(1)=line1
	lint(2)=line2
	planar=.false.
	go to 30
10	continue
30	return
	end
c******************************************************************************
	subroutine setup
c  runtime modification of control parameters
	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 /magxyz/ ef(3),ev(3),tdir(3),azmuth,icompn
	common /pltdev/ jplotr,xboard,yboard,xlimit,ylimit
        character ians*1,ch*1,swt*3
	logical grvswt,magswt
	itty=5
	ierr=0
c
10	print 11
11	format(' setup parameters :'$)
	read(itty,37,end=999) ians
37      format(a)
	call cvc(ians,2)
	n=leftj(ians)
	ch=ians(1:1)
	if(.not.(ch.eq.'d' .or. ch.eq.'m' .or. ch.eq.'p' .or. ch.eq.'s' 
     1  .or. ch.eq.'r')) then
	print *,' setup parameters are :'
	print *,' datum print, magnetic direction'
	print *,' plot parameters, summary, or return'
	if(ierr.ge.2) return
	ierr=ierr+1
	go to 10
	endif
	ierr=0	
c
c  datum print
	if(ians(1:1).eq.'d') then
	if(grvswt) then
	  print 20
20	  format(' set floating gravity datum ? '$)
	  iave=noyes(l)
	endif
	if(magswt) then
	  print 25
25	  format(' set floating magnetic datum ? '$)
	  jave=noyes(l)
	endif
	go to 10
	endif
c
c  magnetic parameters
	if(ians(1:1).eq.'m') then
	call initm(1)
	go to 10
	endif
c
c  plot parameters
	if(ians(1:1).eq.'p') then
	  print *,' current plot device =',jplotr
	  print *,' do you want to change device ?'
	  if(noyes(ldum).eq.1) then
	    print *,' enter new device #: 5=HP Plotter,8=CGA,9=EGA,10=VGA'
	    read(5,*,end=10,err=10) iplotr
	    if(iplotr.ne.jplotr) then
	      xlimit=999.
	      ylimit=999.
	      jplotr=iplotr
	    endif
	    if(jplotr.ne.1) then
	      print *,' enter 0 for an overall plot size '
	      print *,' or    1 for separate scaling factors '
	      read(5,*,end=10,err=10) ifunc
	      if(ifunc.eq.0) then
	        print *,' enter maximum x,y plot size in inches '
	        read(5,*,end=10,err=10) xlimit,ylimit
	      else
	        print *,' enter x,g,h scaling factors'
	        read(5,*,end=10,err=10) xscale,gscale,hscale
	      endif
	    endif
	  endif
c          if(iplotr.eq.1) then
c            xscale=0.0
c            zscale=0.0
c            gscale=0.0
c            hscale=0.0
c          endif
	print *,' (re)set vertical exaggeration and windows ?'
	if(noyes(ldum).eq.1) then
	  print *,' current vertical exaggeration =',vex
	  print *,' enter new vertical exaggeration for the model plot'
	  read(5,*,end=10,err=10) vex
	  print *,' enter 0 0 to allow windows to default'
	  print *,' model vertical window (min,max) :'
	  read(5,*,end=10,err=10) zzz(1),zzz(2)
	  if (grvswt) print *,' gravity window (min,max)  :'
	  if (grvswt) read(5,*,end=10,err=10) ggg(1),ggg(2)
	  if (magswt) print *,' magnetic window (min,max) :'
	  if (magswt) read(5,*,end=10,err=10) hhh(1),hhh(2)
	endif
	print *,' continue with secondary switches ?'
	if(noyes(ldum).eq.0) go to 10
	  print 32
32	  format(' plot station locations ? '$)
	  istatn=noyes(l)
	  print 33
33	  format(' plot body numbers ? '$)
	  numbod=noyes(l)
	  print 34
34	  format(' plot vertex numbers ? '$)
	  labelv=noyes(l)
	  print 35
35	  format(' enable partial derivative plots'/,
     1 ' in the inversion function ? '$)
	  iparsh=noyes(l)
	go to 10
39	print *,' input error'
	go to 10
	endif
c
	if(ians(1:1).eq.'s') then
	print *
	print *,' datum prints are'
	if(grvswt) then
	if(iave.eq.1) then
	  print *,' floating gravity'
	else
	  print *,' nonfloating gravity'
	endif
	endif
	if(magswt) then
	if(jave.eq.1) then
	  print *,' floating magnetics'
	else
	  print *,' nonfloating magnetics'
	endif
	print *
	print *,' magnetic field parameters are'
	print 40,ef
40	format(' intensity=',1pe15.8,/,' inclination=',0pf7.2,/,
     1 ' declination=',f7.2)
	print *,' profile azimuth=',azmuth
	if(icompn.eq.0) then
	  print *,' total field is calculated'
	else
	  dir='t'
	  if(icompn.eq.1) dir='x'
	  if(icompn.eq.2) dir='y'
	  if(icompn.eq.3) dir='z'
	  print 42,dir
42	  format(1x,a1,' component is calculated')
	endif
	endif
	print *
	print *,'plot parameters are'
	print *,' vertical exaggeration =',vex
	swt='off'
	if(istatn.eq.1) swt='on'
	print *,' station location ',swt
	swt='off'
	if(numbod.eq.1) swt='on'
	print *,' body numbers ',swt
	swt='off'
	if(labelv.eq.1) swt='on'
	print *,' vertice labels ',swt
	swt='off'
	if(iparsh.eq.1) swt='on'
	print *,' partial derivative plots ',swt
	go to 10
	endif
c
	if(ians(1:1).eq.'r') return
	go to 10
999	return
	end
c******************************************************************************
	subroutine initw(xxx,grvswt,magswt)
        common /field/  nf(10),f(5120)
	dimension xxx(2)
	logical grvswt,magswt
	xxx(1)=1.e10
	xxx(2)=-xxx(1)
	if(grvswt) then
	  xxx(1)=f(nf(1))
	  xxx(2)=f(nf(2)-1)
	endif
	if(magswt) then
	  if(f(nf(6)  ).lt.xxx(1)) xxx(1)=f(nf(6))
	  if(f(nf(7)-1).lt.xxx(2)) xxx(2)=f(nf(7)-1)
	endif
	return
	end
c******************************************************************************
	subroutine initnf(nobs,mobs,nf,max,ierr)
	dimension nf(10)
	ierr=0
	if(5*(nobs+mobs).gt.max) go to 3
	nf(1)=1
	do 1 i=2,6
1	nf(i)=nf(i-1)+nobs
	do 2 i=7,10
2	nf(i)=nf(i-1)+mobs
	return
3	print 4,max
4	format(' data requires more than the ',i4,' words allocated')
	ierr=1
	return
	end
c******************************************************************************
	subroutine initl(plenth,pl,nbody,ireset)
	dimension pl(2,50)
	if(plenth.lt.0.0) go to 15
	do 10 j=1,nbody
	if(ireset.eq.1) go to 5
	if(plenth.eq.0.0) go to 10
	if(pl(1,j).ne.0. .or. pl(2,j).ne.0.) go to 10
5	pl(1,j)=plenth
	pl(2,j)=-plenth
10	continue
	return
15	print 20
20	format(' initl: no change, plenth must be >= 0.')
	return
	end
c******************************************************************************
	subroutine initp
	common /pltdev/ jplotr,xboard,yboard,xlimit,ylimit
	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
	character*56 title1,title2,title3
	character*20 titlx,titlz,titlg,titlh
	iplotr=9 
	jplotr=iplotr
	xlimit=999.
	ylimit=999.
	do 1 i=1,2
	xxx(i)=0.0
	zzz(i)=0.0
	ggg(i)=0.0
1	hhh(i)=0.0
	sizea =0.08
	sizet =0.1
	sizel =0.15
	xscale=0.0
	zscale=0.0
	gscale=0.0
	hscale=0.0
	adelx =0.0
	adelz =0.0
	adelg =0.0
	adelh =0.0
	lintx=5
	lintz=5
	lintg=5
	linth=5
	titlx='kilometer'
	titlz='depth km'
	titlg='mGal'
	titlh='nTesla'
	return
	end
c******************************************************************************
	subroutine initm(iset)
c  Magnetic field parameters are stored: intensity, inclination, 
c  declination where inc is deg positive down, and dec is deg positive 
c  east of north.  Unit vector components are: x along profile, z 
c  positive down, and y right cartesian.
	common /magxyz/ ef(3),ev(3),tdir(3),azmuth,icompn
	character ians*1,icom*1
	data d2r/1.745329e-2/
	ireset=iset
	if(ef(1).eq.0.0) go to 11
	if(ireset.eq.0) go to 100
c
5	print 10
10	format(' change earth`s field vector ? '$)
	if(noyes(l).eq.0) go to 20
11	print *,' present earth`s field vector '
	print *,'   intensity=',ef(1)
	print *,' inclination=',ef(2)
	print *,' declination=',ef(3)
	print *,' enter 3 new values '
	read(5,*,end=5,err=5) ef
c
20	print 25
25	format(' change the profile azimuth ? '$)
	if(noyes(l).eq.0) go to 30
	print *,' enter new azimuth'
	read(5,*,end=20,err=20) azmuth
c
30	icompn=0
	print  35
35	format(' do you want total field or a component :'$)
	read(5,37,end=999) ians
37      format(a)
	call cvc(ians,2)
	if(ians(1:1).eq.'t') go to 100
	if(ians(1:1).ne.'c') go to 30
50	print  60
60	format(' x,y, or z component :'$)
	read(5,37,end=30) icom
	call cvc(icom,2)
	if(icom(1:1).eq.'x') icompn=1
	if(icom(1:1).eq.'y') icompn=2
	if(icom(1:1).eq.'z') icompn=3
	if(icompn.eq.0) go to 50	
	print  70,icom(1:1)
70	format(' the calculated field values are now the ',a1, 
     1  ' component',/,
     2 ' referenced to the profile direction as the x axis')
c
100	ei=ef(2)*d2r
	ed=(ef(3)-azmuth)*d2r
	eh=cos(ei)
	ev(1)=cos(ed)*eh
	ev(2)=sin(ed)*eh
	ev(3)=sin(ei)
	do 110 i=1,3
110	tdir(i)=ev(i)
	if(icompn.eq.0) return
	do 120 i=1,3
120	tdir(i)=0.
	if(icompn.eq.1) tdir(1)=1.
	if(icompn.eq.2) tdir(2)=1.
	if(icompn.eq.3) tdir(3)=1.
999	return
	end
