c******************************************************************************
	subroutine curvmn(zg[huge],iqd[huge],b,nc,nr,epsmx,nim,
     1  eps1,dn1,ni)
c   applies minimum curvature equations to the first
c  nc*nr elements of array zg.
c   array iqd contains nc*nr elements which indicate
c  for each mesh location the quadrant where a data
c  value is located. an iqd value of zero indicates
c  no data and -1 locks the present mesh value.
c   array b should contain 6*nc*nr elements used for
c  weighting when iqd is 1 to 4, in the case where
c  iqd is only 0 or -1, b can be of length one.
c   the over-relaxation parameter w increases
c  as the system converges until 1.7 is reached.
c  M. Webring USGS
        integer*2 iqd
	dimension zg(1),iqd(1),b(1)
	data nimn/5/,lmtc/1/
	if(nc.lt.5 .or. nr.lt.5) return
	ni=0
	dn=1.e20
	w=1.3
	eps=0.
	eps1=0.
	epsm=abs(epsmx)
111	continue
	if(ni.ge.nim) go to 72
	eps=0.
c first row
	if(iqd(1))2,1,1
1	zg(1)=(( (2.*(zg(2)+zg(nc+1))-zg(nc+nc+1)-zg(3))*.5 )-
     1 zg(1))*w+zg(1)
2	j1=nc+2
	j2=j1+nc
	if(iqd(2))4,3,3
3	zg(2)=(( (4.*(zg(3)+zg(j1))+2.*zg(1)-zg(4)-zg(j1-1)-
     1 zg(j1+1)-zg(j2))*.16666667 )-zg(2))*w+zg(2)
4	do 6 i=3,nc-2
	j1=i+nc
	j2=j1+nc
	if(iqd(i))6,5,5
5	zg(i)=(( (4.*(zg(i-1)+zg(j1)+zg(i+1))-zg(j1+1)-zg(j1-1)-
     1 zg(j2)-zg(i+2)-zg(i-2))*.14285714 )-zg(i))*w+zg(i)
6	continue
	if(iqd(nc-1))8,7,7
7	i=nc-1
	j1=i+nc
	zg(i)=(( (4.*(zg(i-1)+zg(j1))+2.*zg(i+1)-zg(i-2)-
     1 zg(j1+1)-zg(j1-1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
8	if(iqd(nc))10,9,9
9	j1=nc+nc
	zg(nc)=(( (2.*(zg(j1)+zg(nc-1))-zg(nc-2)-zg(j1+nc))*.5 )-
     1 zg(nc))*w+zg(nc)
c second row
10	if(iqd(nc+1))12,11,11
11	i=nc+1
	j1=i+nc
	zg(i)=(( (4.*(zg(j1)+zg(i+1))+2.*zg(1)-zg(2)-
     1 zg(i+2)-zg(j1+1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
12	if(iqd(nc+2))14,13,13
13	i=nc+2
	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i+1))+4.*(zg(jm)+zg(i-1))-
     1 2.*zg(j1+1)-zg(jm+1)-zg(j1-1)-zg(i+2)-zg(j1+nc))*
     1 5.5555556e-2 )-zg(i))*w+zg(i)
14	do 16 i=nc+3,nc+nc-2
	j1=i+nc
	jm=i-nc
	if(iqd(i))16,15,15
15	zg(i)=(( (8.*(zg(i-1)+zg(j1)+zg(i+1))+4.*(zg(jm))-
     1 2.*(zg(j1-1)+zg(j1+1))-zg(jm-1)-zg(jm+1)-
     1 zg(j1+nc)-zg(i+2)-zg(i-2))*5.263158e-2 )-zg(i))*w+zg(i)
16	continue
	i=nc+nc-1
	if(iqd(i))18,17,17
17	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i-1))+4.*(zg(jm)+zg(i+1))-2.*zg(j1-1)-
     1 zg(jm-1)-zg(j1+1)-zg(i-2)-zg(j1+nc))*5.5555556e-2 )-
     1 zg(i))*w+zg(i)
18	i=nc+nc
	if(iqd(i))20,19,19
19	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(j1)+zg(i-1))+2.*zg(jm)-zg(jm-1)-
     1 zg(i-2)-zg(j1-1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
c rows 3 to nr-2
20	do 39 j=3,nr-2
	i=(j-1)*nc+1
	if(iqd(i))22,21,21
21	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(i+1)+zg(j1)+zg(jm))-zg(j1+nc)-zg(j1+1)-zg(i+2)-
     1 zg(jm+1)-zg(jm-nc))*.14285714 )-zg(i))*w+zg(i)
22	i=i+1
	if(iqd(i))24,23,23
23	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i+1)+zg(jm))+4.*zg(i-1)
     1 -2.*(zg(j1+1)+zg(jm+1))-zg(j1-1)-zg(j1+nc)-zg(i+2)-
     1 zg(jm-nc)-zg(jm-1))*5.2631578e-2 )-zg(i))*w+zg(i)
24	do 35 j2=3,nc-2
	i=i+1
	if(iqd(i))35,25,25
25	j1=i+nc
	jm=i-nc
	d=zg(i)
	if(iqd(i))26,26,27
26	d=(( (8.*(zg(i+1)+zg(i-1)+zg(jm)+zg(j1))-2.*(zg(j1+1)+zg(jm+1)+
     1 zg(jm-1)+zg(j1-1))-zg(j1+nc)-zg(jm-nc)-zg(i-2)-zg(i+2))*
     1 .05 )-d)*w+d
	go to 33
27	ndx=(i-1)*6+1
	b1=b(ndx)
	b2=b(ndx+1)
	b3=b(ndx+2)
	b4=b(ndx+3)
	b5=b(ndx+4)
	b6=b(ndx+5)
	go to (28,29,30,31)iqd(i)
28	bu=b1*zg(jm+1)+b2*zg(jm)+b3*zg(i-1)+b4*zg(j1-1)
	go to 32
29	bu=b1*zg(jm-1)+b2*zg(jm)+b3*zg(i+1)+b4*zg(j1+1)
	go to 32
30	bu=b1*zg(j1-1)+b2*zg(j1)+b3*zg(i+1)+b4*zg(jm+1)
	go to 32
31	bu=b1*zg(j1+1)+b2*zg(j1)+b3*zg(i-1)+b4*zg(jm-1)
32	t=.25*(zg(j1+nc)+zg(i-2)+zg(jm-nc)+zg(i+2))
     1 +.5*(zg(j1-1)+zg(jm-1)+zg(jm+1)+zg(j1+1))-
     1 (zg(j1)+zg(i-1)+zg(jm)+zg(i+1))
	d=(( (bu+b5-t)*b6 )-d)*w+d
33	epsln=d-zg(i)
	if(abs(epsln).lt.abs(eps)) go to 34
	eps=epsln
	ieps=i
34	zg(i)=d
35	continue
	i=i+1
	if(iqd(i))37,36,36
36	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(j1)+zg(i-1)+zg(jm))+4.*zg(i+1)-2.*(zg(j1-1)+
     1 zg(jm-1))-zg(jm+1)-zg(jm-nc)-zg(i-2)-
     1 zg(j1+nc)-zg(j1+1))*5.2631578e-2 )-zg(i))*w+zg(i)
37	i=i+1
	if(iqd(i))39,38,38
38	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(j1)+zg(i-1)+zg(jm))-zg(jm-nc)-zg(jm-1)-zg(i-2)-
     1 zg(j1-1)-zg(j1+nc))*.14285714 )-zg(i))*w+zg(i)
39	continue
c row nr-1
40	i=(nr-2)*nc+1
	if(iqd(i))42,41,41
41	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(jm)+zg(i+1))+2.*zg(j1)-zg(jm-nc)-zg(jm+1)-
     1 zg(i+2)-zg(j1+1))*.16666667 )-zg(i))*w+zg(i)
42	i=i+1
	if(iqd(i))44,43,43
43	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i+1)+zg(jm))+4.*(zg(i-1)+zg(j1))-
     1 2.*zg(jm+1)-zg(jm-1)-zg(jm-nc)-zg(i+2)-
     1 zg(j1+1))*5.5555556e-2 )-zg(i))*w+zg(i)
44	do 46 j=3,nc-2
	i=i+1
	if(iqd(i))46,45,45
45	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i-1)+zg(jm)+zg(i+1))+4.*zg(j1)-
     1 2.*(zg(jm-1)+zg(jm+1))-zg(j1-1)-zg(i-2)-
     1 zg(jm-nc)-zg(i+2)-zg(j1+1))*5.2631578e-2 )-zg(i))*w+zg(i)
46	continue
	i=(nr-1)*nc-1
	if(iqd(i))48,47,47
47	j1=i+nc
	jm=i-nc
	zg(i)=(( (8.*(zg(i-1)+zg(jm))+4.*(zg(j1)+zg(i+1))-2.*zg(jm-1)-
     1 zg(j1-1)-zg(i-2)-zg(jm-nc)-zg(jm+1))*5.5555556e-2 )-
     1 zg(i))*w+zg(i)
48	i=i+1
	if(iqd(i))50,49,49
49	j1=i+nc
	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(jm))+2.*zg(j1)-zg(jm-nc)-zg(jm-1)-
     1 zg(i-2)-zg(j1-1))*.16666667 )-zg(i))*w+zg(i)
c last row
50	i=i+1
	if(iqd(i))52,51,51
51	jm=i-nc
	zg(i)=(( (2.*(zg(i+1)+zg(jm))-zg(i+2)-zg(jm-nc))*.5 )-
     1 zg(i))*w+zg(i)
52	i=i+1
	if(iqd(i))54,53,53
53	jm=i-nc
	zg(i)=(( (4.*(zg(i+1)+zg(jm))+2.*zg(i-1)-zg(i+2)-zg(jm+1)-
     1 zg(jm-nc)-zg(jm-1))*.16666667 )-zg(i))*w+zg(i)
54	do 56 j=3,nc-2
	i=i+1
	if(iqd(i))56,55,55
55	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(i+1)+zg(jm))-zg(i-2)-zg(jm-1)-
     1 zg(jm-nc)-zg(jm+1)-zg(i+2))*.14285714 )-zg(i))*w+zg(i)
56	continue
	i=i+1
	if(iqd(i))58,57,57
57	jm=i-nc
	zg(i)=(( (4.*(zg(i-1)+zg(jm))+2.*zg(i+1)-zg(i-2)-
     1 zg(jm-1)-zg(jm-nc)-zg(jm+1))*.16666667 )-zg(i))*w+zg(i)
58	i=i+1
	if(iqd(i))60,59,59
59	jm=i-nc
	zg(i)=(( (2.*(zg(i-1)+zg(jm))-zg(i-2)-zg(jm-nc))*.5 )-
     1 zg(i))*w+zg(i)
60	if(ni)70,70,71
70	eps1=abs(eps/w)
71	ni=ni+1
	if(eps.eq.0) go to 72
	dn1=abs(eps/w)
	if(dn1.le.epsm .and. ni.ge.nimn) go to 72
	dlam=dn1/dn
	dn=dn1
	if(dlam.gt.1.) go to 74
	if(dlam.lt..8) go to 75
	if(w.ge.1.6) go to 75
	w=w+.1
	go to 75
74	if(iconv.eq.lmtc) go to 76
	iconv=iconv+1
	go to 75
76	w=w-.1*aint(dlam*10.-9.11)
	iconv=0
	if(w.lt.1.)w=1.
75	continue
	go to 111
72	return
	end
c******************************************************************************
	subroutine rowio(n,z[huge],iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end
c ====================================================================
	subroutine subtr
c ====================================================================
c  This subroutine subtracts second input grid from the first 
	character*50 infil1,infil2,outfile
	character id*56,pgm*8
	common/term/iw,ir,dval
	common/specs/xo,dx,yo,dy,nc,nr,nz
	dimension a(2000),b(2000)
	ispec=0
c
1	write(iw,800)
800	format(' Enter original data file name:')
	read(ir,801)infil1
801	format(a50)
	open(10,file=infil1,status='old',form='unformatted',
     &mode='read',err=5)
	call grdchk(10,ispec,1)
	go to 10
5	write(iw,802)
802	format(' Error opening file. Try again')
	go to 1
10	write(iw,803)
803	format(' Enter file to subtract from original:')
	read(ir,801) infil2
	open(11,file=infil2,form='unformatted',status='old',
     &mode='read',err=15)
	call grdchk(11,ispec,1)
	go to 20
15	write(iw,802)
	go to 10
20	write(iw,804)
804	format(' Enter output file name')
	read(ir,801) outfile
	write(iw,805)
805	format(' Enter title')
	read(ir,806) id
806	format(a56)
	pgm=' VARMAG '
	open(13,file=outfile,status='unknown',form='unformatted')
	write(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
c
c  Begin subtracting
c
	do 100 j=1,nr
	call rowio(nc,a,-1,10,13,iend)
	call rowio(nc,b,-1,11,13,iend)
	do 50 i=1,nc
	if(a(i).ge.1.e+30.or.b(i).ge.1.e+30) go to 40
	a(i)=a(i)-b(i)
	go to 50
40	a(i)=dval
50	continue
	call rowio(nc,a,0,0,13,iend)
100	continue
	close(10)
	close(11)
	close(13)
	return
	end
c ======================================================================
	subroutine trim
c ======================================================================
c
c  this subroutine will find rows and columns in a standard gridded file
c  that are entirely dvals, delete them and change row and column count and
c  x0 and y0 if necessary.  Intended to omit flagged borders, so may cause
c  problems if internal rows or columns are totally flagged
c
c       reduce to fit PC; put x and idel in work
c      dimension x(5000),idel(5000)
        dimension x(2000),idel(2000)
        common /work/x,idel,idum(56000)
	character*50 infile,outfile,temp1,temp2
	character id*56,pgm*8
	common/term/iw,ir,dval
      temp1='slice1.tmp'
      temp2='slice2.tmp'
802   format(a50)
5     write(iw,801)
801   format(' Enter name of file to be trimmed')
      read(ir,802) infile
      open(10,file=infile,status='old',form='unformatted',mode='read',
     &err=5)
      read(10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
10    write(ir,803)
803   format(' Enter name of output trimmed file')
      read (ir,802) outfile
      open(11,file=temp1,form='unformatted',status='unknown')
      write(11) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
c
c  row slice
c
      iopt=1
      itest=0
      irow=0
      do 25 i=1,nrow
      call rowio(ncol,x,-1,10,11,iend)
      do 15 j=1,ncol
      if(x(j).le.1.e+30) go to 20
15    continue
      if(itest.ne.0) go to 25
      y0=y0+dy
      go to 25
20    itest=1
      irow=irow+1
      call rowio(ncol,x,0,10,11,iend)
25    continue
      ncol2=ncol
      nrow2=irow
      y02=y0
      x02=x0
      go to 100
c
c  column slice
c
c
50    iopt=2
      close(10,status='delete')
      open(10,file=temp2,status='old',form='unformatted')
      read(10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
c       reduce to fit PC
        do 52 i=1,2000
c      do 52 i=1,5000
52    idel(i)=0
      icol=0
      itest=0
      do 70 i=1,nrow
      call rowio(ncol,x,-1,10,11,iend)
      do 55 j=1,ncol
      if(x(j).ge.1.e+30.and.idel(j).eq.0) go to 55
      idel(j)=1
55    continue
70    continue
      do 75 i=1,ncol
      if(idel(i).eq.0) go to 73
      icol=icol+1
      itest=1
      go to 75
73    if(itest.ne.0) go to 75
      x0=x0+dx
75    continue
      ncol2=icol
      nrow2=nrow
      x02=x0
      y02=y0
      rewind 10
      open(11,file=outfile,form='unformatted',status='unknown')
      go to 105
c
c  read in temporary file and write out with new cols & rows
c
100   close(10)
      close(11)
      open(10,file=temp1,form='unformatted',status='old')
      open(11,file=temp2,form='unformatted',status='unknown')
105   read(10) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      pgm=' sliced '
      write(11) id,pgm,ncol2,nrow2,nz,x02,dx,y02,dy
      do 40 j=1,nrow2
      call rowio(ncol,x,-1,10,11,iend)
      go to (39,33),iopt
33    kk=0
      do 35 k=1,ncol
      if(idel(k).eq.0) go to 35
      kk=kk+1
      x(kk)=x(k)
35    continue
      if(kk.ne.ncol2) print*, 'yipe'
39    call rowio(ncol2,x,0,10,11,iend)
40    continue
      close(11)
      go to (50,200),iopt
200   close(10,status='delete')
300   return
      end
c =======================================================================
	subroutine vmsub
c ======================================================================
c
c  This subroutine of varmag calculates the magnetization of high relief
c  topography from the synthetic terrain effects that will give a residual
c  with minimum correlation with terrain. It operates on the premise that
c  target anomalies are unrelated to topography, and therefore have minimum
c  correlation with such.
c
c  command file contains the following parameters.  Whatever parameters
c  that are missing in the command file will be asked for at runtime.
c  A COMMAND FILE IS OPTIONAL.  All filenames in the command file must
c  be surrounded by single quotes.
c
c  mfile    name of original anomaly grid
c  mtfile   name of synthetic terrain effects grid
c  dfile    name of grid of damping factors.
c  j0file   (optional) name of grid containing variable initial
c	    magnetizations (emu/cc).  Not normally recommended.
c  jxfile   (optional) name of grid defining areas of magnetization
c  	    that are already knwon and therefore should remain
c	    fixed.
c  xJ0 	    constant initial magnetization in emu/cc used to calculate
c	    the syntehtic terrain effects (mtfile).  It is essential 
c	    that this indeed was the magnetization used to calculate
c	    the data in mtfile.  It is not necessary to specify xJ0
c	    if j0file is used instead.
c  nwind    window length on one side in number of grid points.  Must be
c	    an odd number not greater than 21. 
c  thresh   correlation threshold below which the damped correlation
c	    between the residual and syntehtic terrain effects is
c	    assumed negligible.
c  xJmin    (optional) minimum magnetization allowed in output magnetization
c	    grid.  If a calculated magnetization is below xJmin, the
c	    xJmin is assigned to the center point of the window.  Useful if
c	    it is known that no rocks in the area arer strongly reversely
c	    magnetized.
c  jfile    name of output magnetizaton grid.
c
c  Note:  NO DVALS ALLOWED IN INPUT GRIDS
c
c       store arrays in common to save space
        common/work/ om,tm,r,rf,xJo,xJJ,fj,idum(13000)
      dimension om(1000,21),tm(1000,21),r(1000),rf(1000),xJo(1000),
     &xJJ(1000),fj(1000)
      character*50 mtfile,mfile,j0file,jxfile,jfile,dfile,rfile
      character cfile*50,id*56,pgm*8,strg1*50,strg2*50
      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
      pgm='VARMAG'                   
1     if(mfile.eq.' ') then
5       write(iw,801)
800	format(a50)
801     format(' Enter name of original aeromag anomaly grid')
        read(ir,800) mfile
      endif
      open(10,file=mfile,status='old',form='unformatted',mode='read',
     &err=5)
      call grdchk(10,ispec,1)
      if(mtfile.eq.' ') then
8       write(iw,881)
881	format(' Enter name of synthetic terrain effects file')
        read(ir,800) mtfile
      endif                                       
      open(11,file=mtfile,form='unformatted',status='old',mode='read',
     &err=8)
      call grdchk(11,ispec,1)
      if(j0file.ne.' ') then
	rewind 11
	call idchk(11,j0file,1)
      endif                                       
      open(13,file=dfile,status='old',form='unformatted',mode='read')
      call grdchk(13,ispec,3)
      write(iw,803)
803   format(' Enter output variable-magnetization file name')
      read(ir,800) jfile
      write(iw,804)
804   format(' Enter title for output magnetization file (car ret for
     & default)')
      read(ir,802) id
802   format(a56)
      if(id.eq.' ') then
        call name(mfile,strg1)
        call name(mtfile,strg2)      
        write(id,805) strg1,strg2,thresh,nwind
805	format('*var mag frm ',a10,'&',a10,',thresh=',f3.2,',nwind=',
     &i2)
      endif
      open(16,file=jfile,status='unknown',form='unformatted')
      write(16) id,pgm,nc,nr,nz,xo,dx,yo,dy
      if(iJx.eq.0) go to 25
      open(18,file=jxfile,status='old',form='unformatted',mode='read')
      read(18) id,pgm,nc,nr,nz,xo,dx,yo,dy
25    open(17,file=rfile,status='unknown',form='unformatted')
      call name(jfile,strg1)
      write(id,806) strg1
806   format('*damped correl coefs assoc. with ',a23)
      write(17) id,pgm,nc,nr,nz,xo,dx,yo,dy
      if(inJ0.eq.2) then
        open(14,file=j0file,status='old',form='unformatted',mode='read')
        call grdchk(14,ispec,3)
      endif
c
c  Begin with first window.  Set up window parms first.
c
      write(ir,890)
890   format(' Thinking...'/)
      nbeg=(nwind-0.999999)/2.+1.
      nrend=nr-nbeg+1
      ncend=nc-nbeg+1
      winds=1.e+0/float(nwind*nwind)
c  read data from all but last row of one window
      do 50 j=1,nwind-1
      read(10)dum,(om(i,j),i=1,nc)
50    read(11)dum,(tm(i,j),i=1,nc)
c  set first half of window to dvals (if not inputting jxfile)
      do 55 i=1,nc
      r(i)=dval
      xJJ(i)=dval
55    continue
      if (iJx.ne.0) go to 56
      do 550 i=1,nc
      fj(i)=dval
550   continue
56    if(inJ0.eq.1) then
        do 57 i=1,nc
        xJo(i)=xJ0
57      continue
      endif
      do 60 j=1,nbeg-1
      read(13) dum,(rf(i),i=1,nc)
      go to (560,555),iJx+1
555   read(18) dum,(fj(i),i=1,nc)
      do 557 i=1,nbeg-1
      xJJ(i)=fj(i)
      kn=nc-i+1
      xJJ(kn)=fj(kn)
557   continue
560   write(16)dum,(fj(i),i=1,nc)
      write(17) dum,(r(i),i=1,nc)
      go to (60,59),inJ0
59    read(14) dum,(xJo(i),i=1,nc)
60    continue
c
c  Begin doing windows for whole grid
c
      jr=nwind-1
      do 500 n=nbeg,nrend
      ictr=nbeg
      jr=jr+1   
      if(jr.gt.nwind) jr=1
c  read next row of data grids in window
      read(10)dum,(om(i,jr),i=1,nc)
      read(11)dum,(tm(i,jr),i=1,nc)
      read(13) dum,(rf(i),i=1,nc)
      go to (62,61),iJx+1
61    read(18) dum,(fj(i),i=1,nc)
      do 661 i=1,nbeg-1
      xJJ(i)=fj(i)
      kn=nc-i+1
      xJJ(kn)=fj(kn)
661   continue
62    go to (65,63),inJ0
63    read(14) dum,(xJo(i),i=1,nc)
c
c  start doing all windows along first set of rows 
c
c  set up default values of correlation array and magnetization
65    ic1=ictr-nbeg+1
      ic2=ictr+nbeg-1
      xJ=xJo(ictr)
      if(iJx.eq.0) go to 100
      if(fj(ictr).ge.1.e+38) go to 100
      xJ=fj(ictr)
      r(ictr)=dval
      go to 400
c  find averages of window
100   sumo=0.0
      sums=0.0
      sumt=0.0
      do 200 j=1,nwind
      do 200 i=ic1,ic2
      if(om(i,j).ge.1.e+30.or.tm(i,j).ge.1.e+30) then
      	print*,'dvals not allowed'
      	stop
      endif
      s=om(i,j)-tm(i,j)
      sumo=sumo+om(i,j)
      sums=sums+s
      sumt=sumt+tm(i,j)
200   continue
      aveo=sumo*winds
      aves=sums*winds
      avet=sumt*winds   
c  find variances needed and slope of regression
      sumot=0.0
      sums=0.0
      sumt=0.0
      sumts=0.0
      do 225 j=1,nwind     
      do 225 i=ic1,ic2
      ss=om(i,j)-tm(i,j)-aves
      tt=tm(i,j)-avet
      oo=om(i,j)-aveo
      sumot=sumot+oo*tt
      sums=sums+ss*ss
      sumt=sumt+tt*tt
      sumts=sumts+ss*tt
225   continue  
c
c  Calculate correl coefs for center of window & multiply by damp. factor
c
      denom=sqrt(sumt*sums)
      if(denom.eq.0.0) go to 400
240   r(ictr)=sumts*rf(ictr)/denom
c
c  Check if minimum correl.  If not, adjust magnetization appropriately.
c
c        write (*,*)' line 701 denom,sumt,sums',denom,sumt,sums
      if(abs(r(ictr)).le.thresh) go to 400
c calculate J in terms of ratio J/Ja by linear regression
c   do not allow xJ to go above xJmax or below xJmin
c        write (*,*) ' sumot,sumt',sumot,sumt
      ratio=sumot/sumt
      xJ=xJo(ictr)*ratio
      if(xJ.gt.xJmax) xJ=xJmax
      if(xJ.lt.xJmin) xJ=xJmin
c
c  output magnetization in xJJ array
c
400   xJJ(ictr)=xJ
      ictr=ictr+1
      if(ictr.le.ncend) go to 65
      write(16) dum,(xJJ(k),k=1,nc)
      write(17) dum,(r(k),k=1,nc)
500   continue
c write last completely dvaled rows (unless inputting jxfile)
      do 510 i=1,nc
      r(i)=dval
510   continue
      if(iJx.ne.0) go to 515
      do 512 i=1,nc
      fj(i)=dval
512   continue
515   do 520 j=1,nbeg-1
      if(iJx.ne.0) read(18) dum,(fj(i),i=1,nc)
      write(16) dum,(fj(i),i=1,nc)
      write(17) dum,(r(i),i=1,nc)
520   continue
      close(16)
      close(17)
      close(10)
      close(11)
      close(13)
      close(14)
      close(18)         
      return
      end
      subroutine namemc(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
        subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c       assigns values to proper variables
c       variables are passed to program varmag through common blocks
c       numr=position in the array var where real variables start
c       numa=position in the array var where arrays start
c       nnvar=number of variables in program varmag
c
        parameter (nnvar=16,numr=4)
        character*6 pvar,var(nnvar)
        character*56 tvar,kvar,cfmt
        logical chv
        character*56 title
        character*50 mtfile,mfile,j0file,jxfile,dfile,cfile,ifile,
     1  rfile,jfile
        character*50 vfile,ofile
        common/vmparms/xJ0,nwind,thresh,xJmin,xJmax
        common/vmfile/mtfile,mfile,j0file,jxfile,jfile,dfile,rfile,
     1  cfile
        common/parms/ vfile,ifile,ofile,narea,title,iopt
        data var/'nwind','narea','iopt',
     1  'xJ0','thresh','xJmin','xJmax',
     2  'vfile','ifile','ofile','title',
     3  'mtfile','mfile','j0file','jxfile','dfile'/
c
        numa = 17
c
c     inum is set to one only to prevent a compilier
c     warning about a variable not being used.
c     it must be removed if any arrays would be used
c     in the future.
c
      inum=1
        do 190 i = 1,nnvar
        if (pvar .ne. var(i)) goto 190
        if (.not. chv) then
c
c       non character value
c       right justify the number in kvar
c
        m = 57 - nn
        im = m - 1
        kvar(m:56) = tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write (cfmt,50) im,nn
50      format ('(',i2,'x,i',i2,')')
        read (kvar,cfmt) jvar
        else
c
c       real
c
        write (cfmt,60) im,nn
60      format ('(',i2,'x,g',i2,'.0)')
        read (kvar,cfmt) xvar
        endif
        endif
        goto (101,102,103,104,105,106,107,108,109,110,
     1  111,112,113,114,115,116),i
101     nwind = jvar
        goto 200
102     narea = jvar
        goto 200
103     iopt = jvar
        goto 200
104     xJ0 = xvar
        goto 200
105     thresh = xvar
        goto 200
106     xJmin = xvar
        goto 200
107     xJmax=xvar
        goto 200
108     vfile = tvar(1:nn)
        goto 200
109     ifile = tvar(1:nn)
        goto 200
110     ofile = tvar(1:nn)
        goto 200
111     title = tvar(1:nn)
        goto 200
112     mtfile = tvar(1:nn)
        goto 200
113     mfile = tvar(1:nn)
        goto 200
114     j0file = tvar(1:nn)
        goto 200
115     jxfile = tvar(1:nn)
        goto 200
116     dfile = tvar(1:nn)
        goto 200
190     continue
        write (*,*) ' error in namelist - ',pvar,
     1  'variable not included'
        stop
200     nvar = i
        return
        end

