c  CHESS
c  Combination of programs CHESSBOARD coded by Lin Cordell, PRECHESS
c  coded by Tien Grauch and MCHESS coded by Tom Hildenbrand, with a later
c  modification by Cordell.
c  Program uses the chessboard method of Lin Cordell, modified with
c  Hildenbrand's sliding-rule filter.
c  Reference to the method:
c     Cordell, Lindrith, 1985, Techniques, applications, and 
c        problems of analytical continuation of New Mexico
c        aeromagnetic data between arbirary surfaces of very
c        high relief [abs.]:  Proceedings of the International
c        Meeting on Potential Fields in Rugged Topography,
c        Institute of Geophysics, University of Lausanne,
c        Switzerland, Bulletin No. 7, p. 96-99.	 
c
c  The chessboard method interpolates between parallel levels or
c  surfaces in order to estimate level-to-surface, surface-to-level or
c  surface-to-surface continuation.  (Because it is an approximate method,
c  users desiring surface-to-surface continuation that involves parallel
c  surfaces will achieve the same but faster results if 
c  program NFFTFIL is used in up/downward continuation mode.  The
c  amount to continue is the distance between the parallel surfaces.)
c
c  The level-to-surface option of the chessboard method is straightforward,
c  but the other two options require fooling the program.  In these
c  two cases, the program constructs a surface grid p to enter for the
c  surface-to-level operation in the following way:
c                    p = 2*l - surf,
c  where l=level to continue to, surf=surface where data are now.
c   OR for the surface-to-surface operation in the following way:
c	             p = newsurf - oldsurf + l,
c  where the data are being continued from oldsurf to newsurf, and l
c  is a level somewhere in between the surfaces to fool CHESS with.
c
c  The sliding rule filter is a taperd cutoff which affects progressively
c  lower wave numbers (longer wave lenghts) with increasing downward
c  continuation.  The fft spectrum is multiblied by a filter weight having
c  a value of 1 at wavelengths >= w1 = 1.5*z + (n+1)*del, and tapers linearly
c  to 0 at wavelength w2 = z + n*del, where del is the grid interval; z is the
c  depth of downward continuation, and n is 1,2 or 3, depending on whether
c  "light", "moderate", or "heavy" filtering is selected.  All spectral values
c  above w2 are set to zero.  No filtering is applied with upward continuation
c  or if "no filtering" is selected.
c
c  Some judgement is called for in using this program, as with too much
c  filtering data continued downward can be made to appear smoother that
c  data continued upward!  In critical cases or where the range of elevation
c  is larger that the grid interval, it is advisable to supply program
c  CHESS with prefiltered coefficients, using, for example, program F_STRIP,
c  with treatment of wrap-around provided for by program PREP.  Also, if
C  magnetic data are to eventually be converted to pseudo-gravity,
c  the pseudo-gravity coefficients can be supplied to provide more stable down-
c  -ward continuation.  To supply your own coefficients, obtain a fftfil.cof
c  file from program FFTFIL, and select the "presupplied coefficients"
c  option.  Note that the continuation is done in the frequency domain, using
c  parts of program FFTFIL, and nomenclature follows that of program FFTFIL.

c
      dimension work(2048),a1(2,1024,16),fz(1024,25)
      dimension z(1024)
c      character crfile*12,iunit*1
      character*50 tname,coname,fname,flnam,fl2nam
      common/bmain[near]/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
      common/parm1[near]/fname,coname,zlev(25),nlev,inocal,barlev,ifac
        common/parm2/tname
      data dvaltst/1.0e+28/,ndim/1024/,dval/0.1701412E+39/
	ir=5
	iw=6
	kr=11
c sm is slop allowed when comparing grid specs
	sm=.000001
	print*                              
5 	write(6,888)
888	format(' 1 - level to surface continuation'/
     1 ' 2 - surface to level continuation'/' 3 - ',
     2 'surface to surface continuation'/5x,'Enter option: ',$)
	read*,iopt
	if(iopt.le.0.or.iopt.gt.3) go to 5
	go to (7,8,8), iopt
7     write(iw,803)
803	format(' Enter input magnetic data (observed on level): '$)
	read(ir,10) fname
	open(17,file=fname,status='old',form='unformatted',mode ='read')
      read(17)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
	close(17)
	if(ncol.gt.ndim) then
		write(iw,810) ndim
810		format(' Max no. cols = ',i4)
		stop
	endif
      write(iw,804)
804	format(' Enter surface grid that you wish to continue to: '$)
      read(5,10) flnam
10    format(a50)
	open(10,file=flnam,status='old',form='unformatted',mode='read')
      read(10)id,pgm,ncol2,nrow2,nz,xo2,dx2,yo2,dy2
	if(abs(xo2-xo).gt.sm.or.abs(yo2-yo).gt.sm) go to 995
	if(ncol2.ne.ncol.or.nrow2.ne.nrow) go to 995
	if(abs(dx2-dx).gt.sm.or.abs(dy2-dy).gt.sm) then
	  write(iw,807) flnam
807   format(' WARNING--dx and dy of file may not match input mag:'
     1 /a50)
	endif
	write(iw,811)
811     format(' Enter observ. level (in same units as surface ',
     1 'elev data): ',$)
	read*,barlev
      go to 110
8     write(iw,823)
823	format(' Enter input magnetic data (observed on surface): '$)
	read(ir,10) fname
	open(17,file=fname,status='old',form='unformatted',mode='read')
      read(17)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
	close(17)
	if(ncol.gt.ndim) then
		write(iw,810) ndim
		stop
	endif
      write(iw,805)
805	format(' Enter surface grid that data are observed on: ',$)
      read(5,10) flnam
	open(10,file=flnam,status='old',form='unformatted',mode='read')
      read(10)id,pgm,ncol2,nrow2,nz,xo2,dx2,yo2,dy2
	if(abs(xo2-xo).gt.sm.or.abs(yo2-yo).gt.sm) go to 995
	if(ncol2.ne.ncol.or.nrow2.ne.nrow) go to 995
	if(abs(dx2-dx).gt.sm.or.abs(dy2-dy).gt.sm) then
	  write(iw,807) flnam
	endif
	if(iopt.eq.3) go to 110
	write(iw,812)
812	format(' Enter level to continue to (in same units as observ.',
     1 ' surf. data): ',$)
	read*,barlev
110	write(6,801) 
801	format(' Enter the constant to convert surf. data into',
     1 ' the same units as dx & dy '/' with z positive up',
     1 ' (e.g. .001 converts meters a.s.l. to km) : ',$)
	read*,conv
	call minmax(10,ncol,nrow,z,conv,xmin,xmax,xmean)
	write(6,891) xmin,xmax,xmean
891	format(' Surface minimum= ',g14.6,' maximum= ',g14.6,
     1 ' mean= ',g14.6/
     2 '  in same units as dx & dy'/)
	rewind 10
      read(10)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
	go to (15,15,12),iopt
12	write(iw,806)
806	format(' Enter file of surface to continue to: ',$)
       read(5,10) fl2nam
	open(11,file=fl2nam,status='old',form='unformatted',mode='read')
       read(11)id,pgm,ncol2,nrow2,nz,xo2,dx2,yo2,dy2
	if(abs(xo2-xo).gt.sm.or.abs(yo2-yo).gt.sm) go to 995
	if(ncol2.ne.ncol.or.nrow2.ne.nrow) go to 995
	if(abs(dx2-dx).gt.sm.or.abs(dy2-dy).gt.sm) then
	  write(iw,807) fl2nam
	endif
	write(6,801) 
	read*,conv3
	call minmax(11,ncol2,nrow2,work,conv3,xmin2,xmax2,xmean2)
	write(6,891) xmin2,xmax2,xmean2
c set arbitrary barlev between the two surfaces (part of fooling schess)
	txmin=amin1(xmin,xmin2)
	txmax=amax1(xmax,xmax2)
	barlev= txmin + (txmax - txmin) * 0.5e0
	rewind 11
        read(11)id,pgm,ncol2,nrow2,nz,xo2,dx2,yo2,dy2
15	barlev=barlev*conv
	barlv2=barlev*2.e0
	write(iw,817)
817	format(' Enter no. of levels for interpolation: ',$)
	read(ir,*) nlev
	if(nlev.lt.7.or.nlev.gt.25) then
		write(iw,818)
818		format(' Max no. levels=25, minimum=7')
		go to 15
	endif
	write(iw,815)
815	format(' How should data be filtered in the downward ',
     1  'continuation?'/,
     2  5x,' 0 - no filtering',/5x,' 1 - slight lowpass',/5x,
     3  ' 2 - moderate lowpass',/5x,' 3 - severe lowpass'/,
     4  ' Enter filtering factor: ',$)
	read(ir,*) ifac
	ifac=ifac-1
	write(iw,816)
816	format(' no. of grid points by which to augment rows '
     1 '& cols in FFT (nadd): ',$)
	read(ir,*) nadd
	write(iw,813)
813	format(' Enter output file name: '$)
	read 10,tname
	write(iw,814)
814	format(' Enter title: '$)
	read 802,id
802	format(14a4)
	open(12,file='mchessz.tmp',status='unknown',form='unformatted')
       write(12)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
	go to (20,30,130),iopt
c
c  LEVEL TO DRAPE
c
20	do 25 j=1,nrow
	read(10) dum,(z(i),i=1,ncol)
	do 22 i=1,ncol
	if(z(i).ge.dvaltst) go to 22
	z(i)=z(i)*conv
22      continue
	write(12) dum,(z(i),i=1,ncol)
25	continue
	go to 300
c
c  DRAPE TO LEVEL
c
30	xmin=dvaltst
	xmax=-dvaltst
	do 100 j=1,nrow
	read(10) dum,(z(i),i=1,ncol)
	do 50 i=1,ncol
	if(z(i).ge.dvaltst) go to 50  
	z(i)=barlv2-(z(i)*conv)
	if (z(i).lt.xmin) xmin=z(i)
	if (z(i).gt.xmax) xmax=z(i)
50      continue
	write(12) dum,(z(i),i=1,ncol)
100	continue
	go to 300
c
c  DRAPE TO DRAPE
c
130	xmin=dvaltst
	xmax=-dvaltst
c find xmin, xmax of new surface while creating it
	do 200 j=1,nrow
	read(10) dum,(z(i),i=1,ncol)
	read(11) dum,(work(i),i=1,ncol)
	do 150 i=1,ncol
	if(z(i).ge.dvaltst.or.work(i).ge.dvaltst) then
		z(i)=dval
		go to 150
	endif
	z(i)=(work(i)*conv3)-(z(i)*conv)+barlev
	if(z(i).lt.xmin) xmin=z(i)
	if(z(i).gt.xmax) xmax=z(i)
150      continue
	write(12) dum,(z(i),i=1,ncol)
200	continue
	go to 300
995	print*,'Grids mismatch'
300	close(11)
	close(10)
      close(12)
        call schess(ndim,xmin,xmax,z,work,a1,fz)
      stop
      end
c ***************************************************************
      subroutine minmax(iunit,ncol,nrow,z,conv,zmin,zmax,zmean)
c  finds min and max of file in converted units and stores converted
c  values back into z
c  the file is already attached to iunit and the header has been read
      dimension z(ncol)
      zmax=-1.e38
      zmin=1.e38
c
c  calculate min, max and mean and convert at same time
c
      sumz=0.0
      en=0.0
      do 20 j=1,nrow
      read(iunit,end=80,err=70) dum,z
      do 10 i=1,ncol
c calculate mean
      if(z(i).ge.1.e+38) go to 10
      z(i)=z(i)*conv
      sumz=sumz+z(i)
      if(z(i).gt.zmax) zmax=z(i)
      if(z(i).lt.zmin) zmin=z(i)
      en=en+1.0
10    continue
20    continue
	if(en.eq.0.0) stop 'found no valid points'
      zmean=sumz/en
	return
70 	print*,'Error in reading the file'
	stop
80	print*,'Found end of file, check pass of ncol & nrow'
	print*,'Has the header been read?'
	stop
	end
c****************************************************************
c    automatic version of chessboard of lin cordell
c    programmed by tom hildenbrand
        subroutine schess(ndim,xmin,xmax,z,work,a1[huge],fz[huge])
      dimension a1(2,ndim,16),work(2*ndim),title(14),pgmd(2)
      dimension fz(ndim,25),z(ndim)
      character*50 tname,coname,fname
	character crfile*12
      common/bmain[near]/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
      common/parm1[near]/fname,coname,zlev(25),nlev,inocal,barlev,ifac
      common/parm2/ tname
      data lnri/16/,dval/1.0e+28/
   20 format(a50)
      if(barlev.lt.xmin) go to 134
      zlev(2)=xmin-0.001
      cmax=barlev
      if(barlev.lt.xmax)cmax=xmax
      xcal=float(nlev-3)
      xlevl=cmax-xmin+0.01
      xdel=xlevl/xcal
      zlev(1)=zlev(2)-xdel/3.
      inocal=nlev-1
      do 138 i=3,nlev-2
      zlev(i)=zlev(i-1)+xdel
138   continue  
      go to 139
134   cmax=xmax
      zlev(2)=(xmin-barlev)/2.+barlev
      zlev(1)=barlev
      xcal=float(nlev-4)
      xlevl=xmax-xmin+0.0001
      xdel=xlevl/xcal
      inocal=1
      zlev(3)=xmin
      do 135 i=4,nlev-2
      zlev(i)=zlev(i-1)+xdel
135   continue                       
139   if(barlev.lt.xmax) go to 137
      zlev(nlev)=xdel+barlev
      zlev(nlev-1)=barlev
      go to 140
137   dif=xmax-barlev
      if(dif.gt.xdel/2.) go to 149
      n=nlev-2
      zlev(nlev-1)=barlev+xdel
      zlev(nlev)=barlev +2.*xdel
      inocal=n
      go to 146
149   zlev(nlev)=xmax+xdel
      zlev(nlev-1)=xmax+0.001
      if(inocal.eq.1) go to 140
      n=nlev-1
144   n=n-1
      inocal=n
      dif=abs(barlev-zlev(n))
      if(dif.lt.xdel/2.or.n.eq.3) go to 146
      go to 144
146   zlev(n)=barlev
140   write(iw,151) (zlev(j),j=1,nlev)
151   format(' levels are as follows:'/,5(5(g15.4,1x),/))
c  read header record of input file.
  180 open(11,status='old',form='unformatted',file=fname)
c  note that x & y have been swicthed from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
      read(kr)title,pgm,n2,n1,nz,yop,dy,xop,dx
      no=inocal+9
      write(crfile,161)no
161   format('mchess',i2,'.tmp')
      open(21,status='unknown',form='unformatted',file=crfile)
      write(21)title,pgm,n2,n1,nz,yop,dy,xop,dx
      do 160 i=1,n1
      read(kr)dum,(z(j),j=1,n2)
      write(21)dum,(z(j),j=1,n2)
160   continue
      close(11)
      close(21)
      open(11,status='old',form='unformatted',file=fname)
      read(kr)title,pgm,n2,n1,nz,yop,dy,xop,dx
      pgm(1)='auto'
      pgm(2)='ches'
      nx=n1
      ny=n2
      l=lnri+1
      lnri21=lnri/2+1
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nx+2*nadd.
c  m=no. of rows, l=no. from 9-16, k=interger.
      n1=n1+2*nadd
  190 l=l-1
      if(l.lt.lnri21) go to 260
      mr=n1/l+0.0000001
      k=1
      idiv=2
  200 iquot=mr/idiv+0.0000001
      if(iquot.lt.idiv) go to 210
      k=k+1
      mr=iquot
      go to 200
  210 k=k+1
      m=l*2**k
  220 mtest=l*2**(k-1)
      if(mtest.lt.n1) go to 230
      k=k-1
      m=mtest
      if(k.eq.0) go to 230
      go to 220
  230 lnxa=m-n1
      if(l.ne.lnri) go to 250
      nxa16=lnxa
  240 nri=l
      nxa=lnxa
      go to 190
  250 if(lnxa.ge.nxa) go to 190
      go to 240
  260 n1=n1+nxa
c  check to see if row block size of 16 will be more efficient
      n116=n1-nxa+nxa16
      ntest=0.9*n116
      if(ntest.gt.n1.or.n116.gt.1028) go to 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 n2=n2+2*nadd
      if(n1.gt.1028.or.n2.gt.1028) go to 320
      nxa=nxa+2*nadd
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
c  input parameters for fft are set: call main program.
	call cfftfil(a1,work,nx,nri,nxa,n1,n2,id2)
c  Interpolates a function between levels, given a sequence of grids at 
c  successive vertical levels. Z direction + downward, input grids
c  in monotonic-increasing-Z sequence only.  All grid specs identical.
c  Extra top and bottom level grids are created.
      open(10,status='old',form='unformatted',file='mchessz.tmp')
      read(10) title,pgmd,ncol,nrow,nz,x0,dx,y0,dy
      do 4 i=1,nlev
      no=9+i
      write(crfile,888)no
888   format('mchess',i2,'.tmp')
      iunit=i+11
c      print*,iunit,no
      open(unit=iunit,status='old',form='unformatted',file=crfile)
      read(iunit) title,pgmd,ncol1,nrow1,nz1,x01,dx1,y01,dy1
4     continue
      open(11,status='unknown',form='unformatted',file=tname)
      write(11) id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      do 400 irow=1,nrow
      read(10)dum,(z(k),k=1,ncol)
      do 19 n=1,nlev
      iunit=n+11
19      read(iunit) dum,(fz(j,n),j=1,ncol)
      do 18 j=1,ncol
      if(z(j).gt.dval) go to 18
      if(fz(j,4).lt.dval) go to 665
      z(j)=fz(j,4)
      go to 18
665   zj=z(j)
c  Note: Assume 1<iclose<nlev, no checks !
        
      do 7 i=2,nlev-1
      if(zlev(i).gt.zj) go to 8
      if(zlev(i).lt.zj) go to 7
      z(j)=fz(j,i)
      go to 18
7     continue
8     iclose=i-1
        IF (ICLOSE.LT.2.OR.ICLOSE.GE.NLEV) THEN
        WRITE(*,*)' CAUTION: addressing undefined level values.'
        WRITE(*,*)'          results will be invalid'
        WRITE(*,*)'          in SCHESS, iclose = ',ICLOSE
C       Added warning--see note above
c       With valid data, this should never occur.  MRM 10-13-88
        ENDIF
      Z2=zlev(iclose)
      x1=zlev(iclose-1)-Z2
      x3=zlev(iclose+1)-Z2
      x4=zlev(iclose+2)-Z2
      x1sq=x1*x1
      x3sq=x3*x3
      x4sq=x4*x4
      f2=fz(j,iclose)
      if(x1.eq.0..or.x3.eq.0..or.x4.eq.0.) go to 620
      q1=(fz(j,iclose-1)-f2)/x1
      q3=(fz(j,iclose+1)-f2)/x3
      q4=(fz(j,iclose+2)-f2)/x4                
      denom=(x3*x4sq-x4*x3sq)-x1*(x4sq-x3sq)+x1sq*(x4-x3)
      if(denom.eq.0.) go to 630
      a=(q1*(x3*x4sq-x4*x3sq)-x1*(q3*x4sq-q4*x3sq)
     1+x1sq*(q3*x4-q4*x3))/denom
      b=(q3*x4sq-q4*x3sq-q1*(x4sq-x3sq)+x1sq*(q4-q3))/denom
      c=(q4*x3-q3*x4-x1*(q4-q3)+q1*(x4-x3))/denom
           x=zj-Z2
           xsq=x*x
      z(j)=f2+a*x+b*xsq+c*xsq*x
      go to 18
620     write(6,640)irow,j
640     format(' x1,x3,or x4=0.: row & col= ',2i4)
      z(j)=0.
      go to 18
630   write(6,650)irow,j
650   format(' denom=0.: row & col= ',2i4)
      z(j)=0.
18    continue
      write(11)dum,(z(k),k=1,ncol)
400    continue
900   continue
      close(10,status='delete')
      close(11)
      do 21 i=1,nlev
      iunit=i+11
21    close(unit=iunit,status='delete')
      go to 350
c  errors resulting in job abortion.
  280 write(iw,290)
  290 format(' #parameter error detected---case run aborted')
      go to 350
  320 write(iw,330)nx,ny,n1,n2
  330 format(' #no. of extended rows or columns exceeds 1028:'/
     1' input no. of rows and columns='2i4/
     2' no. of rows and columns required for filtering=',2i4,/)
  340 close(11)
  350 write(iw,360)
  360 format(' end of job')
      stop
      end
