c =================================================================
c  PROGRAM VARMAG
c =================================================================
c  This program employs the variable-magnetization terrain-correction
c  method that is based on the premise that magnetic sources of
c  interest are often geomtrically unrelated to terrain.  The method
c  finds the magnetization that gives a magnetic-field residual with
c  minimum correation to terrain effects for a window of data within
c  a grid of magnetic-field values.  By repeating the calculation
c  for windows covering the entire grid, a grid of variable magneti-
c  zation values is produced which is combined with topography to 
c  calculate a magnetic-terrain correction.
c
c  The procedure used on one window of data is as follows:
c
c  1. Choose an initial magnetization (xJ0) for the study area.
c
c  2. Calculate the magnetic effects of topogrpahy using this xJ0 
c     (accomplished with R. Blakely's program PFMAG3D, USGS
c     Open-File Report 81-298).
c
c  3. Calculate a residual (s) by subtracting the results of step
c     #2 (xJ0*t) from the original aeromagnetic data (f)
c
c  4. Calculate the correlation coefficient (r) between the residual
c     (s) and results of step #2 (xJ0*t).
c
c  5. Multiply the correlaton coefficient (r) by the damping factor
c     (rf) for that grid point (see explanation under subroutine
c     dfactr).
c
c  6. Is the damped correlation coefficient below a user-given
c     correlation threshold?  If so, assign xJ0 to the center grid point
c     of the window, and go to next window (step #3).
c
c  7. If not, calculate xJ by a least-square fit between f and xJ0*t
c     and assign xJ to the center grid point.  Proceed with next window
c     at step #3.
c
c  VARMAG can be run interactively or with a command file.  It is 
c  organized in terms of functions; each function has specific parameters,
c  input, and output associated with it.  Following is a list and brief
c  description of available functions.  Only the first two letters of the
c  function name are recognized by the program.
c
c	dfactrs - calculates damping factors.
c	correction - runs terrain-correction procedure
c	plug - plugs flagged areeas of missing data using minimum
c	  curvature to extrapolate.
c	trim - trims off flagged boundary from magnetization grids.
c	output - saves output grids of damping factors and/or of
c	  damped correlation coefficients
c	input - inputs file of saved damping factors to the program.
c	residual - subtracts theoretical terrain effects from original
c	  anomaly data. (The theoretical terrain effects must be
c	  calculated separately using PFMAG3D.)
c	save - save all known parameters for terrain correction in a
c	  new command file.
c	change - change terrain-correction parameter values.
c	edit - edit or create magnetization grids for input.
c	jxfile - input areras of magnetization that remain fixed during
c	  terrain correction (areas where magnetization is fairly well
c	  known, for instance).
c	mean - remove mean from input grid.  Recommended before terrain
c	  correction
c	type - type values of parameters.
c	list - quick list of functions available.
c	help - brief description of functiions and where to get help.
c
c  Internal switches are as follows:
c
c   inJ0  0 - not set
c	  1 - input constant initial magnetization
c	  2 - input grid of initial magnetizations
c
c   ispec 0 - not set
c	  1 - a grid has been read into program; can now check
c		other grids against it.
c
c   idf   0 - not set
c	  1 - damping factors have been calculated
c	  2 - damping factors were input from a saved file
c
c   jout  0 - a jfile has not been calculated
c	  1 - a jfile has been calculated
c
c   iJx   0 - do not input jxfile
c	  1 - input jxfile
c
c   igrd  0 - plug or trim functions not processed
c         1 - plug or trim function has been processed
c
c  V.J.S. (Tien) Grauch, USGS, February 1986
c
    	character*50 mtfile,mfile,j0file,jxfile,jfile,dfile,rfile
        character cfile*50,func*2,answer*1,strg*20
        character*20 dummy
	common/term/iw,ir,dval
	common/vmparms/xJ0,nwind,thresh,xJmin,xJmax
	common/vmfile/mtfile,mfile,j0file,jxfile,jfile,dfile,rfile,cfile
	common/specs/xo,dx,yo,dy,nc,nr,nz
	common/switch/inJ0,ispec,idf,jout,iJx,igrd
c	namelist/parms/xJ0,nwind,thresh,xJmin,xJmax,mtfile,mfile,
c     &j0file,jxfile,dfile
c	data dval/1.701412E38/,xJ0/1.e+38/,nwind/0/,thresh/-1.0/
c	data xJmin/-1.e+38/,xJmax/1.e+38/,inJ0/0/,ispec/0/,idf/0/
c	data jout/0/,iJx/0/,igrd/0/,iw/6/,ir/5/
        dval=1.701412E38
        xJ0 = 1.0E38
        nwind = 0
        thresh = -1.0
        xJmin=-1.0E38
        xJmax=1.0E38
        inJ0=0
        ispec=0
        idf=0
        jout=0
        iJx=0
        igrd=0
        iw=6
        ir=5

c dval is a large floating point number approx. equal to 10**37
	mtfile=' '
	mfile=' '
	j0file=' '
	jxfile=' '
	jfile=' '
	dfile=' '
	rfile='rfile.tmp'
c
801	format(a2)
802	format(a1)
803	format(a50)
813	format(a20)
c
c Ask for a function.
c
1	mv=0
	write(iw,800)
800	format(/' Function : ',$)
	read(ir,801) func
c ***************************************************************
c              function to exit program
c ***************************************************************
	if(func.eq.'st'.or.func.eq.'ST'.or.func.eq.'ex'.or.func.eq.
     &'EX'.or.func.eq.'qu'.or.func.eq.'QU') go to 100
c ***************************************************************
c        change parameter values - function change
c ***************************************************************
	if(func.eq.'ch'.or.func.eq.'CH') then 
5	  write(iw,805)
805	  format(' 1. Initial mag (xJ0 or j0file)'/
     & ' 2. Window length (nwind)'/' 3. Correl. threshold (thresh)'/
     & ' 4. Min/max magnetiztn allowed (xJmin, xJmax)'/
     & ' 5. Return to function level'/5x,' Enter number of parameter
     & to change : ',$)
	  read(ir,*) np
	  if(np.lt.1.or.np.gt.5) go to 1
	  go to (10,20,25,30,1), np
c ---change initial mag.----
10	  write(iw,806)
806	  format(' Constant initial mag or grid (c or g)?')
	  read(ir,802) answer
	  if(answer.eq.'g') then
	  	inJ0=2
12		write(iw,807)
807		format(' Enter name of init. mag grid (emu/cc)')
	  	read(ir,803) j0file
	  	open(14,file=j0file,status='old',form='unformatted',
     1 mode='read')
c     & mode='read',err=12)
	  	call grdchk(14,ispec,2)
	  else
	  	inJ0=1
	  	write(iw,808)
808		format(' Enter constant init. magnetiztn (emu/cc)')
	  	read(ir,*) xJ0
	  endif
15	  if(mv.ne.0) go to 40
	  write(iw,809)
809	  format(' Enter new synthetic terrain effects file that
     & corresponds to new init. mag.')
	  read(ir,803) mtfile
	  open(12,file=mtfile,status='old',form='unformatted',
     1  mode='read')
c     & mode='read',err=15)
	  call grdchk(12,ispec,1)
	  rewind 12
	  call idchk(12,j0file,2)
	  write(iw,899)
899       format(' Remember to recalculate damping factors if
     &terrain is different')
	  go to 5
c ----change window length----
20	  write(iw,810)
810	  format(' Enter nwind in no. of grid pts (max=21,must be odd)')
	  read(ir,*) nwind
	  if(nwind.le.0.or.nwind.gt.21) go to 20
	  if(mod(float(nwind),2.0).eq.0) go to 20
	  if(mv.ne.0) go to 42
	  go to 5
c ----change correl threshold----
25	  write(iw,811)
811	  format(' Enter correl. threshold (0<thresh<=1)')
	  read(ir,*) thresh
	  if(thresh.le.0.0.or.thresh.gt.1.0) go to 25
	  if(mv.ne.0) go to 43
	  go to 5
c ----change min & max J allowed----
30	  strg=' '
	  write(iw,812)
812	  format(' Enter new min. allowable J (optional)')
	  read(ir,813) strg     
        dummy=strg
        call shiftr(dummy,20)
	  if(strg.ne.' ') read(dummy,'(g20.0)') xJmin
	  write(iw,814)
814	  format(' Enter new max. allowable J (optional)')
	  read(ir,813) strg 
        dummy = strg
        call shiftr(dummy,20)
	  if(strg.ne.' ') read(strg,'(g20.0)') xJmax
	  go to 5
	endif
c **************************************************************
c       calculate damping factors - function dfactor
c **************************************************************
	if(func.eq.'df'.or.func.eq.'DF') then
	  if(idf.eq.0) dfile='dfile.tmp'
	  call dfactr(ispec)
	  idf=1
	  go to 1
	endif
c **************************************************************
c      plug or trim jfile grid - functions plug and trim
c **************************************************************
	if(func.eq.'pl'.or.func.eq.'PL'.or.func.eq.'tr'.or.
     & func.eq.'TR') then
	  if(func.eq.'pl'.or.func.eq.'PL') then
	  	call plug
	  else
	  	call trim
	  endif
	  igrd=1
	  go to 1
	endif
c ***************************************************************
c      terrain correction procedure - function correction
c ***************************************************************
	if(func.ne.'co'.and. func.ne.'CO') goto 8818
32	  write(iw,887)
887	  format(' Enter input command file name (optional)')
          read(ir,803) cfile
          if(cfile.eq.' ') go to 35
          open(9,file=cfile,form='formatted',
     &mode='read',status='old',err=32)
        call namemc(9)
c          read(9,parms)
          close(9)      
35	  if(dfile.eq.' ') then
		write(iw,817)
817		format(' Must calculate or input damping factors
     & first')
	  	go to 1
	  endif
	  if(inJ0.eq.0) then
		if(xJ0.lt.1.e+38.and.j0file.eq.' ') then
		  inJ0=1
		  go to 40
		endif
		if(j0file.ne.' ') then
		  inJ0=2
		  go to 40
		endif	  
	  	mv=1
	  	go to 10
	  endif
40	  if(nwind.eq.0) then
	  	mv=1
	  	go to 20
	  endif
42	  if(thresh.lt.0.0) then
	  	mv=1
	  	go to 25
	  endif
          if(jxfile.ne.' ') iJx=1
43	  call vmsub
	  jout=1
	  igrd=0
	  go to 1
8818    continue
c *****************************************************************
c  save all known parm values in new command file - function save
c *****************************************************************
	if(func.eq.'sa'.or.func.eq.'SA') then
	  answer=''''
	  nu=9
	  write(iw,818)
818	  format(' Enter output command file name')
	  read(ir,803) cfile
	  open(nu,file=cfile,status='unknown',form='formatted')
	  write(nu,897)
897	  format(' $parms')
44	  write(nu,888) answer,mtfile,answer,answer,mfile,answer,
     & nwind,thresh,xJmin,xJmax
888	  format(' mtfile=',a1,a50,a1,','/' mfile=',a1,a50,a1,','/
     & ' nwind=',i2,',thresh=',f5.2,',xJmin=',g14.7,',xJmax=',g14.7)
	  if(inJ0.eq.2) then
		write(nu,889) answer,j0file,answer
889		format(' j0file=',a1,a50,a1,',')
	  endif
	  if(inJ0.eq.1) then
		write(nu,890) xJ0
890		format(' xJ0=',g14.7,',')
	  endif
	  if(iJx.eq.1) write(nu,891) answer,jxfile,answer
891	  format(' jxfile=',a1,a50,a1,',')
	  if(dfile.ne.'dfile.tmp'.and.dfile.ne.' ') write(nu,892)
     & answer,dfile,answer
892	  format(' dfile=',a1,a50,a1,',')
	  if(nu.eq.9) write(nu,894)
894	  format(' $')
	  if(nu.eq.9) close(nu)
	  go to 1
	endif
c *****************************************************************
c       edit or create a magnetization file - function edit
c *****************************************************************
	if(func.eq.'ed'.or.func.eq.'ED') then
	  call edtmag
	  go to 1
	endif
c *****************************************************************
c       type parameter values on terminal - function type
c *****************************************************************
	if(func.eq.'ty'.or.func.eq.'TY') then
	  nu=iw
	  answer=' '
	  go to 44
	endif
c *****************************************************************
c subtract calculated frm observed to residual - function residual
c *****************************************************************
	if(func.eq.'re'.or.func.eq.'RE') then
	  call subtr
	  go to 1
	endif
c *****************************************************************
c   output damping factors or damped corr coefs - function output
c *****************************************************************
	if(func.eq.'ou'.or.func.eq.'OU') then
	  write(iw,821)
821	  format(' Want to output file of damping factors?')
	  read(ir,802) answer
	  if(answer.eq.'y'.or.answer.eq.'Y') then
	  	if(idf.eq.0) then
	  	  write(iw,819)
819		  format(' Need to calc damping factors first')
	  	  go to 1
	  	endif
		write(iw,822)
822		format(' Enter output damp. factors file')
		read(ir,803) cfile
		call copyf(dfile,cfile)
	  	dfile=cfile
	  endif
	  if(jout.eq.0) go to 1
	  write(iw,823)
823	  format(' Want to output damped correl coefs from last
     & correction?')
	  read(ir,802) answer
	  if(answer.eq.'y'.or.answer.eq.'Y') then
	  	write(iw,882)
882		format(' Enter output damp. correl coef file')
		read(ir,803) cfile
		call copyf(rfile,cfile)
	        rfile=cfile
	  endif
	  go to 1
	endif
c *****************************************************************
c    input areas of J that should remain fixed - function jxfile
c *****************************************************************
	if(func.eq.'jx'.or.func.eq.'JX') then
	  if(iJx.eq.0) then
45		iJx=1
	  	write(iw,824)
824		format(' Name of file containing areas of J that
     & should remain fixed (jxfile)')
	  	read(ir,803) jxfile
	  	open(15,file=jxfile,status='old',form='unformatted',
     & mode='read',err=45)
	  	call grdchk(15,ispec,2)
	  	go to 1
	  else
	  	write(iw,826)
826		format(' Discontinue inputting any areas of fixed J?')
	  	read(ir,802) answer
	  	if(answer.ne.'y'.and.answer.ne.'Y') go to 45
	  	iJx=0
	  	go to 1
	  endif
	endif
c *******************************************************************
c              input damping factors - function input
c *******************************************************************
	if(func.eq.'in'.or.func.eq.'IN') then
	  incount=0
36	  write(iw,827)
827	  format(' Enter damping factors file name')
	  read(ir,803) dfile
	  open(13,file=dfile,status='old',form='unformatted',mode='read',
     &err=37)
	  call grdchk(13,ispec,0)
	  idf=2
	  go to 1
37	  incount=incount+1
	  if(incount.ge.2) then
		write(iw,828)
828		format(' I guess you''d better calculate them instead')
		go to 1
	  else
		  write(iw,898)
898		  format(' Try again')
		  go to 36
	  endif
	endif
c *******************************************************************
c               remove mean from grid - function mean
c *******************************************************************
	if(func.eq.'me'.or.func.eq.'ME') then
	  call mean
	  go to 1
	endif
c *******************************************************************
c               type help instructions - function help
c *******************************************************************
	if(func.eq.'he'.or.func.eq.'HE') then
	  write(iw,831)
831	  format(' This program performs a terrain-correction procedure
     & for aeromagnetic data'/' to give a variable magnetization grid.
     &   Also performs related functions.'/' The method is NOT automatic,
     & so PLEASE read the written documentation before'/' continuing!
     &  Following is a brief description of functions (only the first'/
     & ' 2 letters will be recognized).'//
     & ' dfactor - calculates damping factors'/
     & ' correction - runs terrain-correction procedure'/
     & ' plug - plugs flagged areas using minimum curvature'/
     & ' trim - trims off flagged boundaries of output var-mag file'/
     & ' output - saves damping factors or damped correl coefs'/
     & ' input - inputs damping factors file to program'/
     & ' residual - subtracts theor. terrain effects from original
     & anomaly data'/
     & ' save - save all known parameters in new command file'/
     & ' change - change parameter values'/
     & ' edit - edit or create magnetization grids'/
     & ' jxfile - input areas of magnetization that remain fixed
     & during terrain correct.'/
     & ' mean - remove mean from grid'/
     & ' type - type values of parameters'/
     & ' list - quick list of functions available'//
     & ' **NOTE: Theoretical terrain effects should be calculated
     & with pfmag3d**')
	go to 1
	endif
c *******************************************************************
c           write quick list of functions - default
c *******************************************************************
	write(iw,830)
830	format(/15x,'QUICK LIST OF FUNCTIONS (sees only 1st 2 chars)'//
     & ' dfactors',12x,'correction',10x,'plug',16x,'trim'/
     & ' output (files)',6x,'residual',12x,'save (parms)',8x,
     & 'change (parms)'/
     & ' edit (mag)',10x,'jxfile',14x,'type (parms)',8x,'help'/
     & ' list',15x,' input',15x,'mean',12x,'stop, quit, or exit')
	go to 1
c *****************************************************************
c     END OF PROGRAM - delete temporary files where applicable
c *****************************************************************
100     if(idf.eq.1) then
	  open(13,file='dfile.tmp',status='old')
	  close(13,status='delete')
	endif
	if(jout.eq.1) then
	  open(17,file='rfile.tmp',status='old')
	  close(17,status='delete')
	endif
	if(igrd.eq.0.and.jout.eq.1) then
	  write(iw,840) jfile
840	  format('     WARNING: ',a50,/' should be plugged or
     & trimmed before using to calculate'/' theoretical terrain
     & effects.')
	endif
	stop
	end
c ===========================================================================
      subroutine copyf(infile,outfile)
c ===========================================================================
c  This subroutine copies the infile to the outfile (a rename routine)
      dimension x(2000)
c       Put x in the common area, work, to samve space
        common/work/x,idum(58000)
      character*50 infile,outfile
      character id*56,pgm*8
c
      open(10,file=infile,form='unformatted',status='old',mode='read')
      open(11,file=outfile,form='unformatted',status='unknown')
      read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      write(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
      do 50 j=1,nr
      call rowio(nc,x,1,10,11,iend)
50    continue
      close(10)
      close(11)
      return
      end
c ======================================================================
 	subroutine dfactr(ispec)
c ======================================================================
c
c  Subroutine in program varmag that calculates damping factors to be
c  multiplied to correlation coefficients. Routine incorporates the
c  horizontal gradient algorithm from R. Simpson that approximates a
c  derivative at a grid point by dividing the difference between the
c  two neighbors by 2*grid interval.  The damping reduces the magnitude
c  of the correlation in areas where the variation of the data is small.
c  The factor is empirical and for this program relies only on the
c  variation of one file (for varmag it is variation of synthetic
c  terrain effects).  The factor is
c
c	1 - exp(-G),  where G is
c
c                  g
c           -----------------
c              (ave. of g)
c
c  g=magnitude of horiz gradient 
c
c  Tien Grauch  April 1985
c
      character*50 mtfile,dfile,strg,mfile,j0file,jxfile,jfile,rfile
      character*50 cfile
      character id*56,pgm*8
      dimension a(2000),b(2000),c(2000),grad(2000),aa(2000),bb(2000)
      dimension cc(2000)
      common/term/iw,ir,dval
      common/vmfile/mtfile,mfile,j0file,jxfile,jfile,dfile,rfile,cfile
      common/specs/xo,dx,yo,dy,nc,nr,nz
c       Put variables in common area, work, to save space
        common/work/a,b,c,grad,aa,bb,cc,idum(46000)
      data pgm/'dfactor'/
c
1       write(iw,800)
800	format(' Name of input synthetic terrain effects file:')
        read (ir,801)mtfile
801     format(a50)
      open(12,file=mtfile,status='old',form='unformatted',mode='read',
     &err=5)
      call grdchk(12,ispec,1)
      go to 7
5     write(iw,802)
802   format(' file not found or wrong type.  Try again')
      go to 1
c
7     open(13,file='dfact1.tmp',form='unformatted',status='unknown')
      write(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
c
c  Begin calculation of gradient of grid
c      
      call rowio(nc,b,-1,12,13,iend)
      call rowio(nc,c,-1,12,13,iend)
      nc1=0
      sum1=0.0
c
      do 10 i=1,nc
 10   grad(i)=dval
c
      do 200 j=2,nr-1
      do 120 i=1,nc
      a(i)=b(i)
 120  b(i)=c(i)
      call rowio(nc,c,-1,12,13,iend)
      do 40 i=1,nc
 40   grad(i)=dval
      do 100 i=2,nc-1
      if(b(i-1).ge.1.0e+38.or.b(i+1).ge.1.0e+38) go to 100
      if(a(i).ge.1.0e+38.or.c(i).ge.1.0e+38) go to 100
      dzdx=(b(i+1)-b(i-1))/(2.0*dx)
      dzdy=(c(i)-a(i))/(2.0*dy)
      grad1=sqrt(dzdx*dzdx+dzdy*dzdy)
      nc1=nc1+1
      sum1=sum1+grad1
      grad(i)=grad1
 100  continue
c
      call rowio(nc,grad,0,12,13,iend)
 200  continue
c  Fill in last row
      do 510 i=1,nc
 510  grad(i)=dval
      call rowio(nc,grad,0,12,13,iend)
c       apr 15,1988 changed close to unit 12
        close (12)
c      close(10)
c
c  Calculate averages of gradients
c
      ave1=sum1/float(nc1)
c
c  Read multiplied gradients back in, calc. factor
c
      rewind 13
      read(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
      open(20,file=dfile,form='unformatted',status='unknown')
      call name(mtfile,strg)
      write(id,805) strg
805   format('*damping factors for ',a35)
      write(20) id,pgm,nc,nr,nz,xo,dx,yo,dy
c  First row is all dvals
      do 550 i=1,nc
      a(i)=dval
550   continue
      call rowio(nc,a,0,13,20,iend)
      do 700 j=2,nr-1
      b(1)=dval
      b(nc)=dval
      call rowio(nc,grad,-1,13,20,iend)
      do 650 i=2,nc-1
      if(grad(i).ge.1.e+30) then
	b(i)=dval
	go to 650
      endif
      G=grad(i)/ave1
      b(i)=1.0-exp(-G)
650   continue
      call rowio(nc,b,0,13,20,iend)
700   continue
c  Last row is all dvals
      call rowio(nc,a,0,13,20,iend)
      close(13,status='delete')
      close(20)
600   return
      end
c ==========================================================================
	subroutine edtmag
c =========================================================================
c
c  create or edit areas of polygonal boundaries defined in the command file.
c  The areas are given a specific magnetization value as specified in command
c  file.  The magnetization files are used in conjunction with VARMAG and
c  PFMAG3D.        iopt defines whether to create a new jxfile or j0file
c  or whether to edit an old magnetization file.  Only narea, and definition
c  of the areas are required to appear in the command file; vfile must exist.
c  Other parameters will be asked for as needed.  Can save the parameters in a 
c  new command file after answering all the questions.
c     Max 200 vertices, 50 areas.
c
c  Most of the code, including the crucial 'inside' function were coded by
c  Mike Webring, USGS.  Modified for VARMAG by Tien Grauch January 1986.
c
c namelist parameters
c
c vfile	file containing x,y coordinate pairs that define the polygon areas.
c	The coordinates must be in data units (km, ft, etc.--whatever dx
c	and dy are in).
c	vfile is read in free-field format with one coordinate pair per line.
c
c narea	number of polygonal areas defined following the namelist
c
c
c ****************EXAMPLE**************************************************
c
c -------------------Set up------------------------------------------------
c
c 3  1--------2---------------3       Two separate areeas referenced to the
c    |         \         2    |       same coordinate system.  
c 2  |    1     7-----8-------6       
c    |           \              
c 1  4------------5
c
c    1  2  3  4  5  6  7  8  9
c
c -----------------vertice file--------------------------------------------
c
c 1 3
c 4 3
c 9 3
c 1 1        x,y coordinates of the vertices, in the order as labeled
c 5 1            above
c 9 2
c 5 2
c 7 2
c
c ------------------command file-------------------------------------------
c &parms
c narea=2,vfile='vfile.dat',jxfile='jxfile.grd',title='Create a jxfile',
c iopt=1
c &
c 1.2e-03     << magnetization of the area in emu/cc
c 5           << number of vertices in the next line
c 1 2 7 5 4   << list of vertices defining a polygon in clockwise order
c 2.5e-03
c 5
c 2 3 9 8 7
c
c --------------------------------------------------------------------------
c
	common/term/iw,ir,dval
	common/vmparms/xJ0,nwind,thresh,xJmin,xJmax
	common/specs/xo,dx,yo,dy,nc,nr,nz
	dimension list(200),xvert(200),yvert(200),xv(200),yv(200),xJ(50)
	dimension z(2000)
	character*50 ifile,ofile,vfile,cfile,tmp(2),cfile2
	character id*56,title*56,pgm*8,answer*1
	logical inside
	common /parms/ vfile,ifile,ofile,narea,title,iopt
c
	iopt=0
	narea=0
	vfile=' '
	title=' '
	ifile=' '
	ofile=' '
	nask=0
	write(iw,800)
800	format(' enter command filename (car ret to exit):'$)
	read(ir,801) cfile
801	format(a50)
	if(cfile.eq.' ') return
	open(unit=9,file=cfile,status='old',form='formatted',mode='read')
        call namemc(9)
c	read(9,parms)
c  check parms; prompt if some missing
	if(narea.eq.0) stop 'narea must be specified'
	if(narea.gt.50) stop 'max 50 areas allowed'
	if(iopt.eq.0) then
1	  write(iw,802)
802	  format(/' 1 - Create new grid having certain known J areas'/
     &    ' 2 - Create new grid of J0''s'/' 3 - Edit existing J grid'/
     &	  5x,'Enter option : ',$)
	  read(ir,*) iopt
	  nask=nask+1
	  if(iopt.lt.1.or.iopt.gt.3) go to 1
	endif
c
	if(vfile.eq.' ') then
	  nask=nask+1
	  write(iw,803)
803	  format(' Enter name of vertice file')
	  read(ir,801) vfile
	endif
	open(unit=12,file=vfile,form='formatted',status='old',mode='read')
	nvert=0
5	nvert=nvert+1
	if(nvert.gt.200) stop ' vertex list > 200 entries'
	read(12,*,end=6) xvert(nvert),yvert(nvert)
	go to 5
6	nvert=nvert-1
	close(12)
c
	if(iopt.ne.3) then
	  if(dx.eq.0.0.or.nc.eq.0.or.nr.eq.0) then
		nask=nask+1
		write(iw,804)
804		format(' Enter existing grid from which to get proper
     & grid specs')
		read(ir,801) tmp(1)
		open(10,file=tmp(1),form='unformatted',status='old',
     & mode='read')
		read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
		close(10)
	  endif
	else
c
	  if(ifile.eq.' ') then
		nask=nask+1
		write(iw,805)
805		format(' Enter J file to edit')
		read(ir,801) ifile
	  endif
	  open(13,file=ifile,status='old',form='unformatted',mode='read')
	  read(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
	endif
	write(iw,806)xo,yo,dx,dy,nc,nr
806	format(' Grid specs are:'/' xo=',g14.7,' yo=',g14.7,' dx=',
     &  g14.7,' dy=',g14.7,/,' nc=',i4,' nr=',i4)
c
	if(ofile.eq.' ') then
	  nask=nask+1
	  write(iw,807)
807	  format(' Enter output J file name')
	  read(ir,801) ofile
	endif
	tmp(1)='edtmag.1'
	if(narea.eq.1) tmp(1)=ofile
	tmp(2)='edtmag.2'
	pgm='edtmag'
	iswitch=1
	open(15,file=tmp(1),status='unknown',form='unformatted')
c
	if(title.eq.' ') then
	  nask=nask+1
	  write(iw,808)
808	  format(' Enter title for output')
	  read(ir,809) title
809	  format(a56)
	endif
	write(15) title,pgm,nc,nr,nz,xo,dx,yo,dy
c
	if(iopt.eq.2.and.xJ0.ge.dval) then
	  nask=nask+1
	  write(iw,810)
810	  format(' Enter default J0 (A/m) for areas not given by 
     & vfile')
	  read(ir,*) xJ0
	endif
	if(iopt.eq.2) write(iw,811) xJ0
811	format(' Using J0 of ',g14.7,' for unspecified areas')
c save parameters in new command file if desired
	ic=0
	if(nask.ge.2) then
	  write(iw,812)
812	  format(' Want to save these parameters in a new command file?')
	  read(ir,813) answer
813	  format(a1)
	  if(answer.eq.'y'.or.answer.eq.'Y') then
		ic=1
		write(iw,814)
814		format(' Enter new command file name: ',$)
		read(ir,801) cfile2
		open(19,file=cfile2,status='unknown',form='formatted')
		write(19,815) ifile,vfile,ofile,title,narea,iopt
815		format(' $parms'/' ifile=''',a50,''','/
     & ' vfile=''',a50,''','/
     & ' ofile=''',a50,''','/,' title=''',a50,''','/,
     & ' narea=',i2,',iopt=',i1/' $')
	  endif	  
	endif
c read info from bottom of command file and edit values into areas
	k=0
20	k=k+1
	read(9,*) xJ(k)
	read(9,*) ncrnr
	if(ic.eq.1) then
	  write(19,*) xJ(k)
	  write(19,*) ncrnr
	endif
	if(ncrnr.le.2) then
	  close(15,status='delete')
	  close(13)
	  close(9)
	  if(ic.eq.1) close(19)
	  stop 'need at least 3 pts to define an area'
	endif
	read(9,*) (list(i),i=1,ncrnr)
	if(ic.eq.1) write(19,*) (list(i),i=1,ncrnr)
c compile an array of vertices for this area
	do 30 i=1,ncrnr
	xv(i)=xvert(list(i))
30	yv(i)=yvert(list(i))
c guarantee that xmn>=xo,xmx<=x2,etc.  
	xmx=xv(1)
	xmn=xv(1)
	ymx=yv(1)
	ymn=yv(1)
	do 50 i=2,ncrnr
	if(xv(i).gt.xmx) xmx=xv(i)
	if(xv(i).lt.xmn) xmn=xv(i)
	if(yv(i).gt.ymx) ymx=yv(i)
	if(yv(i).lt.ymn) ymn=yv(i)
50	continue
	xmn=amax1(xo,xmn)
	x2=dx*float(nc-1)+xo
	xmx=amin1(x2,xmx)
	ymn=amax1(yo,ymn)
	y2=dy*float(nr-1)+yo
	ymx=amin1(y2,ymx)
c  calculate rows and cols assoc. with xmn,xmx,ymn,ymx
	ix=int((xmn-xo)/dx)
	xo2=xo+dx*float(ix)
	iy=int((ymn-yo)/dy)
	yo2=yo+dy*float(iy)
	nc2=int((xmx-xmn)/dx)+1
	nr2=int((ymx-ymn)/dy)+1
	if(ix.lt.0) stop 111
	if(iy.lt.0) stop 222
	if(nc2.gt.nc) stop 333
	if(nr2.gt.nr) stop 444
	if(nc2.lt.0) nc2=1
	if(nr2.le.0) nr2=1
	ixs=ix+1
	ixe=ix+nc2
	iys=iy+1
	iye=iy+nr2
c
	go to (60,65,110), iopt
c
c  create new magnetization grid (options 1 and 2)
c    option 1
60	constant=dval
	go to 70
c    option 2
65	constant=xJ0
70	do 75 i=1,nc
	z(i)=constant
75	continue
c write default values beginning rows not in the area
	if(iys.gt.1) then
	  do 80 j=1,iys-1
80	  call rowio(nc,z,0,13,15,ie)
	endif
c assign default value if not inside the area
	yp=yo2
	do 90 j=iys,iye
	xp=xo2
	do 85 i=ixs,ixe
	z(i)=constant
	if(inside(ncrnr,xv,yv,xp,yp)) z(i)=xJ(k)
85	xp=xp+dx
	call rowio(nc,z,0,13,15,ie)
90	yp=yp+dy
c write default in ending rows not inside the area
	if(iye.lt.nr) then
	  do 95 i=1,nc
95	  z(i)=constant
	  do 100 j=1,nr-iye
100	  call rowio(nc,z,0,13,15,ie)
	endif
	close(15)
c close if narea=1, set up for next run if not
	if(narea.eq.1) go to 200
	iswitch=2
	if(narea.eq.2) tmp(2)=ofile
	open(13,file=tmp(1),status='old',form='unformatted')
	open(15,file=tmp(2),status='unknown',form='unformatted')
	iopt=3
	read(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
	write(15) id,pgm,nc,nr,nz,xo,dx,yo,dy
	go to 20
c option 3 (includes opt 1 & 2 for k > 1)
c  read & write beginning rows not in the area
110	if(iys.gt.1) then
	  do 120 j=1,iys-1
	  call rowio(nc,z,-1,13,15,ie)
120	  call rowio(nc,z,0,13,15,ie)
	endif
c assign new J value if inside the area
	yp=yo2
	do 130 j=iys,iye
	call rowio(nc,z,-1,13,15,ie)
	xp=xo2
	do 125 i=ixs,ixe
	if(inside(ncrnr,xv,yv,xp,yp)) z(i)=xJ(k)
125	xp=xp+dx
	call rowio(nc,z,0,13,15,ie)
130	yp=yp+dy
c read & write ending rows not inside the area
	if(iye.lt.nr) then
	  do 150 j=1,nr-iye
	  call rowio(nc,z,-1,13,15,ie)
150	  call rowio(nc,z,0,13,15,ie)
	endif
c finalize files or setup for next area
	close(15)
	if(k.eq.1) then
	  close(13)
	else
	  close(13,status='delete')
	endif
	if(k.eq.narea) go to 200
	if(iswitch.eq.1) then
	  if(k.eq.narea-1) tmp(2)=ofile
	  iswitch=2
	  open(13,file=tmp(1),status='old',form='unformatted')
	  open(15,file=tmp(2),status='unknown',form='unformatted')
	else
	  if(k.eq.narea-1) tmp(1)=ofile
	  iswitch=1
	  open(13,file=tmp(2),status='old',form='unformatted')
	  open(15,file=tmp(1),status='unknown',form='unformatted')
	endif
	read(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
	write(15) id,pgm,nc,nr,nz,xo,dx,yo,dy
	go to 20
c end of program
200	close(9)
	if(ic.eq.1) close(19)
	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
	if(d.lt.0.0) inside = .not.inside
5	continue
	return
	end
c ========================================================================
	subroutine grdchk(iunit,ispec,iclose)
c ========================================================================
c This subroutine checks the grid specs of a grid against previously looked
c at grids.  Rejects grid if nz > 1, dx .ne. dy or nc > 2000 also.
c
c  iunit =  fortran unit
c  ispec =  0 if haven't looked at a grid yet,
c	    1 if grid specs are in common
c  iclose = 1 leave file open with header read if grid OK
c	    2 close file upon return if grid OK
c	    3 read header and leave open without checking grid
c
	common/term/iw,ir,dval
	common/specs/xo,dx,yo,dy,nc,nr,nz
	character id*56,pgm*8
c
	read(iunit) id,pgm,nc2,nr2,nz2,xo2,dx2,yo2,dy2
	if(iclose.eq.3) return
	if(nz2.ne.1.or.dx2.ne.dy2) then
		write(iw,800)
800		format(' nz must be 1 and dx must equal dy')
		close(iunit)
		stop
	endif
c set common values if this is first grid checked
	if(ispec.eq.0) then
		if(nc2.gt.2000) then
			write(iw,801)
801			format(' No. of cols greater than 2000')
			close(iunit)
			stop
		endif
		nz=nz2
		xo=xo2
		yo=yo2
		dx=dx2
		dy=dy2
		nc=nc2
		nr=nr2		
		ispec=1
		go to 25
	endif
	slop=0.00001
	if(abs(xo-xo2).gt.slop) go to 50
	if(abs(yo-yo2).gt.slop) go to 50
	if(abs(dx-dx2).gt..001*dx) go to 50
	if(nc.ne.nc2) go to 50
	if(nr.ne.nr2) go to 50
25	if(iclose.eq.1) return
	close(iunit)
	return
50	write(iw,803)
803	format(' Grid specs don''t match those of grids in use')
	close(iunit)
	stop
	end
c ======================================================================
	subroutine idchk(iunit,file,iclose)
c ======================================================================
c This subroutine checks to see if init. mag. is correct in a terrain
c effects file (mtfile)
c
c  iunit=   fortran unit of the file to check
c  file=    name of j0file (iopt=1) or mtfile (iopt=2)
c  iclose=1 leave file open with header read upon return
c        =2 close file upon return
c
c Calls subroutine name
c
	character id*56,pgm*8,strg*50,gname*50,file*50
	common/vmparms/xJ0,nwind,thresh,xJmin,xJmax
	common/switch/inJ0,ispec,idf,jout,iJx,igrd
	common/term/iw,ir,dval
c
	read(iunit) id,pgm,nc,nr,nz,xo,dx,yo,dy
c default titles have a * in the first character
	if(id(1:1).ne.'*') then
		if(iclose.eq.2) close(iunit)
		return
	endif
c
c check to see if init. mag correct
c
	if(inJ0.eq.2) then
		call name(file,gname)
		if(id(43:56).eq.gname(1:14)) go to 100
		write(iw,800) file
800		format(' SERIOUS WARNING!  Synthetic terrain effects
     & grid may not match init. mag. grid',/a50)
		write(iw,801) id
801		format(' Title of terrain effects file is ',a56)
	else  
        call shiftr(strg,50)
		read(strg,'(g50.0)') xJtest
		write(*,*) ' xJtest=',xJtest
		if(abs(xJtest-xJ0).lt.0.0001*xJ0) go to 100
		write(iw,802) xJ0
802		format(' SERIOUS WARNING!  Synthetic terrain effects
     & grid may not have used init mag ',g14.7)
		write(iw,801) id
	endif
100	if(iclose.eq.2) close(iunit)
	return
	end
c ====================================================================
	subroutine mean
c ====================================================================
c  This subroutine removes mean from input grid
	character*50 infile,outfile
	character id*56,pgm*8
	common/term/iw,ir,dval
	dimension a(2000)
c       Put a in common area, work
        common/work/ a,idum(58000)
c
1	write(iw,800)
800	format(' Name of grid from which to remove mean (car ret
     & to exit):')
	read(ir,801)infile
801	format(a50)
	if(infile.eq.' ') return
	open(10,file=infile,status='old',form='unformatted',
     &mode='read',err=5)
	read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
	go to 10
5	write(iw,802)
802	format(' Error opening file. Try again')
	go to 1
10	write(iw,804)
804	format(' Enter output file name')
	read(ir,801) outfile
	pgm='MEAN OUT'
	open(13,file=outfile,status='unknown',form='unformatted')
	write(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
c
c  Find mean
c
	sum=0.0
	icount=0
	do 100 j=1,nr
	call rowio(nc,a,-1,10,13,iend)
	do 50 i=1,nc
	if(a(i).ge.1.e+30) go to 50
	icount=icount+1
	sum=sum+a(i)
50	continue
100	continue
	if(icount.eq.0) then
	  write(iw,806)
806	  format(' No valid points in grid')
	  return
	endif
	ave=sum/float(icount)
	rewind 10
	read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
	do 200 j=1,nr
	call rowio(nc,a,-1,10,13,iend)
	do 150 i=1,nc
	if(a(i).ge.1.e+30) go to 150
	a(i)=a(i)-ave
150	continue
	call rowio(nc,a,0,10,13,iend)
200	continue
	write(iw,807) ave
807	format(' Mean of ',g14.7,' removed.')
	close(10)
	close(13)
	return
	end
c =====================================================================
	subroutine name(file,strg)
c =====================================================================
c This subroutine finds the filename from a VAX pathname
c  file =  input pathname
c
c  strg =  output filename
c
	character*50 file,strg
	istart=1
	iend=len(file)
	n=index(file,':')
	if(n.ne.0) then
		istart=n+1
		m=index(file(istart:iend),']')
		if(m.ne.0) then
			istart=istart+m+1
		endif
	endif
c
	strg=file(istart:iend)
	return
	end
c ======================================================================
	subroutine plug
c ======================================================================
c  program megaplug by M. Webring, modified to be a subroutine.
c  program uses the Briggs (1974) minimum curvature algorithm to fill
c  in flagged areas of grids with extrapolated data.
	character*50 infile,outfile
c        character id*56,pgm*8
	common/term/iw,ir,dval
	data idv,jdv/10,11/
37	format(a50)
1	write(iw,800)
800	format(' Enter name of file to be plugged')
	read(ir,37) infile
	open(idv,file=infile,status='old',form='unformatted',
     &mode='read',err=1)
	write(iw,801)
801	format(' Enter output name for plugged magnetization file:')
	read(ir,37) outfile
	open(jdv,file=outfile,status='unknown',form='unformatted')
	call megautp(dval,ir,idv,jdv)
	close(idv)
	close(jdv)
	return
	end
c******************************************************************************
	subroutine megautp(dval,itty,idv,jdv)
c  plug holes using minimum curvature interpolation
c  gridr and curvmn are from the minc program, usgs open file 81-1224
c  M. Webring USGS 
c       changing dimension of work aea for pc, april, 1988
        integer*2 iqd
        common/work/ zg(40000),iqd(40000)
        dimension wz(1000)
c	dimension /work/ zg(250000),iqd(250000),wz(1000)
	character id*56,p*8,p2*8
c       changing dimension of work area for pc, april 1988
        data nwrk/40000/,p2/'m-c plug'/
c	data nwrk/250000/, p2/'m-c plug'/
	rewind idv
	read(idv) id,p,nc,nr,nz,xo,dx,yo,dy
	nn=nc*nr
	if(nn.gt.nwrk) then
	  write (*, 801) nwrk
801	  format(' no. cols * no. rows is >',i6)
	  return
	endif
1	write(*, 802)
802	format(' no. of min. curvature iterations to use (normally 20):'$)
	read(itty,*,err=1) nim
	do 5 i=1,nn
5	iqd(i)=-1
	ndx=1
	do 20 j=1,nr
	call rowio(nc,zg(ndx),-1,idv,jdv,ie)
	i2=ndx
	do 10 i=1,nc
	if(zg(i2).gt.1.e29) then
	  zg(i2)=dval
	  iqd(i2)=0
	endif
10	i2=i2+1
	if(ie.eq.1) go to 99
20	ndx=ndx+nc
c
	call gridr(nc,nr,zg,wz,dval,ier)
	if(ier.eq.1) then
	  write ( *,'(a)') ' gridr: initialization error'
	  return
	endif
	if(nim.eq.0) go to 30
	eps=0.
	call curvmn(zg,iqd,wz,nc,nr,eps,nim,a,b,ni)
30	write(jdv)id,p2,nc,nr,nz,xo,dx,yo,dy
	ndx=1
	do 50 i=1,nr
	call rowio(nc,zg(ndx),0,idv,jdv,ie)
50	ndx=ndx+nc
99	return
	end
c******************************************************************************
	subroutine gridr(nc,nr,zg[huge],wz,dval,ier)
c  quickly replace dvals with reasonable values.
c  A control net spaced at 'nsep' is generated from ring averages
c  with radii 1, 3, nsep; linear interpolation completes the process.
c  Array 'wz' is at least max(nc,nr).
c  M. Webring USGS
	dimension iw1(3)
	dimension zg(1),wz(1)
	data iw1/3,7,9/
	ier=0
	nn=nc*nr
	nsep=(iw1(3)-1)/2
	ns1=nsep+1
	mxpass=max(nc,nr)/(2*nsep)
	i=1
201	if(i.gt.nn) return
	if(zg(i).ge.dval) go to 200
	i=i+1
	go to 201
c
c  insert control point net
c     distances for ring averages specified by iw1
200	if(nc.lt.nsep+ns1 .or. nr.lt.nsep+ns1) go to 140
	ipass=0
	iflag=1
130	iw=iw1(iflag)
	ihw=(iw-1)/2
	iset=ihw*nc+ihw
131	nass=0
	do 110 jj=ns1,nr-ihw,nsep
	ip=(jj-1)*nc+ns1
	do 100 ii=ns1,nc-ihw,nsep
	if(zg(ip).lt.dval) go to 100
	ips=ip-iset
	it=0
	t=0.0
	do 121 j=1,iw
	ip2=ips
	do 120 i=1,iw
	if(zg(ip2).ge.dval) go to 120
	t=zg(ip2)+t
	it=it+1
120	ip2=ip2+1
121	ips=ips+nc
	if(it-3.lt.0) go to 100
	zg(ip)=t/float(it)
	nass=nass+1
100	ip=ip+nsep
110	continue
	ipass=ipass+1
	if(nass.eq.0 .and. iflag.eq.3) go to 140
	if(ipass.gt.mxpass) go to 140
	if(iflag-2)133,132,131
133	iflag=2
	go to 130
132	iflag=3
	go to 130
140	continue
c
c  fill holes
	inc=nc*nsep
	j=inc+1
	if(nr.lt.ns1) go to 21
	do 20 irow=ns1,nr-nsep,nsep
	call plugm(nc,zg(j),dval)
20	j=j+inc
21	do 28 icol=1,nc
	j=icol
	do 22 k=1,nr
	wz(k)=zg(j)
22	j=j+nc
	call plugm(nr,wz,dval)
	j=icol
	do 24 k=1,nr
	zg(j)=wz(k)
24	j=j+nc
28	continue
c  final check
	do 40 i=1,nn
	if(zg(i).ge.dval) go to 41
40	continue
	return
41	t=0.0
	it=0
	do 42 i=1,nn
	if(zg(i).ge.dval) go to 42
	t=t+zg(i)
	it=it+1
42	continue
	if(it.eq.0) stop ' cannot init grid'
	t=t/float(it)
	write(*,43)t
43	format(' gridr init with',1pe15.5)
	do 44 i=1,nn
	if(zg(i).ge.dval) zg(i)=t
44	continue
	return
	end
c******************************************************************************
	subroutine plugm(n,z[huge],dv)
c  plug holes using linear interpolation
c  M. Webring USGS
	dimension z(n)
	do 1 is=1,n
	if(z(is) .lt. dv) go to 2
1	continue
	return
2	ix=is
3	ix=ix-1
	if(ix.lt.1) go to 4
	z(ix)=z(is)
	go to 3
4	do 5 idv=is,n
	if(z(idv) .ge. dv) go to 6
5	continue
	return
6	is=idv-1
	do 7 ie=idv,n
	if(z(ie) .lt. dv) go to 10
7	continue
	ix=is
9	ix=ix+1
	if(ix.gt.n) return
	z(ix)=z(is)
	go to 9
10	dz=(z(ie)-z(is))/float(ie-is)
	do 11 i=is+1,ie-1
11	z(i)=z(i-1)+dz
	is=ie
	go to 4
	end
c       ************************************************************
        subroutine shiftr(a,ll)
c       added mar,1988 for pc
c       shift list directed real input to right of field for internal read
        character*20 a,c
        character*1 b(20)
        equivalence (c,b)
        c=a
        do 10 i=ll,1,-1
        if (b(i) .ne. ' ') goto 11
10      continue
        i=1
11      ishft=ll-i
        a=' '
        a(ishft+1:ll)=c(1:i)
        return
        end














