c*******************************************************************************
c       program fftfil
c   filter operations employing f.f.t.
c   coded by T.G. Hildenbrand--usgs menlo park 415-329-5303
c   Added function vertot, coded by Hildenbrand and T. Grauch
c   Sept. 20, 1988 (parameters xinc,dec,bdec,binc are
c   required for vertot)
c   (subroutines 'sfftmg' & 'sdummy' coded by R. Watts)
c   structure of command file--
c   Some repairs effected by Hildenbrand Jan 91.
c   More repairs by Lin Cordell 12 Apr 91.
c
c    card 1 : operation in coded name from list below
c              psdmag    pseudo-magnetic transformation
c              psdgrv    pseudo-gravity transformation
c              redpol    reduction of magnetics to the north pole
c              upcont    upward continuation
c              dncont    downward continuation
c              1stver    1st-vertical derivative of the input field
c              2ndver    2nd-vertical derivative of the input field
c              banpas    bandpass filter
c              strike    directional filtering.
c              vertot    change vertical-component data to
c                        total-field
c              nofilt    no filtering
c    card 2 : input file name.ext.
c    card 3 : output file name.ext.
c    card 4 : title (col. 1-56).
c    card 6 : &parms
c    card 6 :  list parameters--example: iopt1=0,iopt2=1,den=.3,....
c    card 7 : &end
c    the following queries are asked after filter operation
c    has been completed. purpose of queries are to continue
c    filtering input file with different filters.
c
c    query 1 : additional filter to be applied? (y or n) - if 'y'
c              the program asks queries 2 thru 6. if 'n' the
c              program starts over by asking for the command segment.
c    query 2 : new operator? format same as in card 1 response.
c    query 3 : new output file name.ext?
c    query 4 : new title?
c    query 5 : parmameter change? (y or n) - if 'n' the filtering
c              computations are started. if 'y' the user enters the
c              parmeters to be changed in a namelist
c              (eg. &parms thet1=-90.,thet2=0.,&)
c              (note the following paramters cannot be changed-nadd,
c              iopt2,idval,ddx,ddy,xo,yo).
c    query 1 is then repeated for additional filter operations.
c   input parameters:
c    iopt1 = 0 no printed output (default: iopt1=0)
c            6 output printed on terminal
c           -1 output printed on disk.
c    iopt2 = -1 no removal of mean from input array (default iopt2=-1)
c             0 remove mean using boundary values
c             1 remove mean and save grid, file name: fftfil.grd.
c    nadd = no. of rows or columns added to each side of grid
c           to reduce the effects of gibbs phenomena (default nadd=0).
c    w1,w2,w3 & w4 - wavelengths used in bandpass filtering.
c           (default w1=w2=0. w3=w4=1.0e+30 infinte wavelength)
c    den - density contrast, gm/cc (default den=1.).
c    bmag - magnetization contrast, gammas (default bmag=1.).
c    dec & xinc - declination and inclination of earth''s field
c          ,in degrees. (default dec=0. xinc=90.)
c    bdec & binc - declination and inclination of magnetization vector.
c                  (default bdec=0. binc=90.)
c    idval = 0 no flagged grid points in input data (default idval=0)
c            1 flagged grid points in input data
c           -1 flagged values removed and locations in file 'flag.loc'.
c    icoef = 1 save fourier coefficients in segment 'fftfil.cof' for
c              later use but perform designated filter operation
c              (note: if data contains flagged values a segment
c              called 'flag.loc' containing their locations is also
c              created in users disk area.)
c            0 fourier coefficients not saved (default icoef=0).
c           -1 fourier coefficients in segment 'fftfil.cof' are used
c              as input. note that if flagged values are present in the
c              data, segment 'flag.loc' is required and 'idval' must be
c              equal to -1. in addition, the parameter 'nadd' must be
c              indentical to its assigned value when the fourier
c              coefficients were saved.
c     z = continuation distance (must be in grid units) (default z=0.)
c         for z>0 downward continuation, for z<0 upward continuation
c    thet1 & thet2 - angles from geographic north that form a pie-slice
c                    filter for directional filtering(-90.ge.thet.le.
c                    +90.;thet2.gt.thet1). (default thet1=0. thet2=90.)
c    istr = -1 reject trends between thet1 and thet2.
c           +1 pass (default istr=1).
c    drive = (RAM) drive letter for temporary files.  By writing temporary
c            files to a RAM drive, calculation speed can be increased by a
c            factor of 2 or more.  See the DOS help file for RAMDRIVE.SYS.
c   (following parameters can be used if spacing and origin
c   of grid want to be changed. note that these changes
c   only occur in the header of the output filtered grid and no
c   subset of the input grid is extracted).
c    ddx - new grid spacing in x-direction.
c    ddy - new grid spacing in y-direction.
c    xo - new origin of rows.
c    yo - new origin of coloumns. (default ddx=ddy=xo=yo=0.)
c
c  Compilation: link with SFFTMG, SFOURT, SFFTFIL
c***********************************************************************
      dimension a1(2,2048,16),work(4096),title(14)
      character*50 tname,coname,fname,conout
      character help*4,blank*2,exit*2,fopr*6,flist(11)*6,blid*4
        character*4 pgm
      character*2 drive
      common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
      common/parm1/fname,coname,iopt1,iopt2,idval,icoef,ddx,ddy
      common/parm2/den,dec,xinc,bmag,bdec,binc,z,thet1,thet2,istr,
     1w1,w2,w3,w4,drive
      equivalence (id(1),blid)
      data icmd/10/,lnri/16/
      data nlist,flist/11,'psdgrv','redpol','upcont','dncont','2ndver',
     1'strike','banpas','psdmag','1stver','vertot','nofilt'/
      data blank/'  '/,exit/'ex'/,help/'help'/
      kr=11
      iw=6
      ir=5
c read command file and header record of input file holding data array.
	print*
   10 write(iw,15)
   15 format(' Enter command file name (car ret to exit): ',$)
      read(ir,20)tname
   20 format(a50)
      if(tname.eq.exit.or.tname.eq.blank) go to 350
c  default values of input parameters.
      iopt1=0
      iopt2=-1
      icoef=0
      nadd=0
      w1=0.
      w2=0.
      w3=1.0e+30
      w4=1.0e+30
      den=1.
      z=1.
      dec=0.0
      xinc=90.
      bmag=1.
      bdec=0.0
      binc=90.
      xo=0.0
      yo=0.0
      ddx=0.0
      ddy=0.0
      thet1=0.0
      thet2=90.
      istr=1
      idval=0
      z=0.0
      drive='  '
c  read command file.
      open(10,status='old',form='formatted',file=tname,share='denywr')
      read(icmd,30)fopr
   30 format(a6)
      if(fopr.ne.blank) go to 50
   40 write(iw,45)
   45 format(' enter operator: ',$)
      read(ir,30)fopr
   50 if(fopr.eq.help) go to 80
      do 60 i=1,nlist
      if(fopr.ne.flist(i)) go to 60
      itype=i
      go to 110
   60 continue
      write(iw,70)
   70 format(' #invalid operator')
   80 write(iw,90)
   90 format(' quick list of options')
      write(iw,100)flist
  100 format(2x,a6)
      go to 40
  110 read(icmd,20)fname
      if(fname.ne.blank) go to 120
      write(iw,115)
  115 format(' enter input file name: ',$)
      read(ir,20)fname
  120 read(icmd,20)coname
      if(coname.ne.blank) go to 130
      write(iw,125)
  125 format(' enter output file name: ',$)
      read(ir,20)coname
  130 read(icmd,140)id
      if(blid.ne.blank) go to 150
      write(iw,135)
  135 format(' enter title: ')
      read(ir,140)id
  140 format(14a4)
  150 call snamel(icmd)
      close(10)
c  test for parameter errors.
      if(z.ge.0.0.and.fopr.eq.'upcont') go to 300
      if(z.le.0.0.and.fopr.eq.'dncont') go to 300
      if(fopr.eq.'nofilt'.and.icoef.eq.0) go to 280
      if(abs(thet1).gt.90..or.abs(thet2).gt.90.) go to 280
      if(thet1.gt.thet2) go to 280
      if(w1.gt.w2.or.w2.gt.w3.or.w3.gt.w4) go to 280
      if(abs(binc).gt.90..or.abs(xinc).gt.90.) go to 280
      if(iopt1)170,180,160
  160 if(iopt1.ne.6) go to 280
      go to 180
  170 iopt1=20
      conout='fftfil.out'
      open(20,status='unknown',form='formatted',file=conout)
c  read header record of input file.
  180 open(11,status='old',form='unformatted',file=fname,share='denywr')
c  note that x & y have been switched 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
      if(icoef.eq.-1.and.fopr.eq.'nofilt')itype=12
      pgm(1)='fftf'
      pgm(2)='il**'
      nx=n1
      ny=n2
      if(ddx.ne.0.0)dx=ddx
      if(ddy.ne.0.0)dy=ddy
      if(xo.ne.0.0)xop=xo
      if(yo.ne.0.0)yop=yo
      xo=xop
      yo=yop
      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.2048) go to 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 n2=n2+2*nadd
      write(iw,275)nri
  275 format(' blocking rows =',i3)
      if(n1.gt.2048.or.n2.gt.2048) 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 sfftfil(a1,work,nx,nri,itype,nxa,n1,n2,id2)
      if(iopt1.eq.20)close(20)
      go to 10
c  errors resulting in job abortion.
  280 write(iw,290)
  290 format(' #parameter error detected---case run aborted')
      go to 10
  300 write(iw,310)
  310 format(' #sign of z does not conform to continuation operator')
      go to 10
  320 write(iw,330)nx,ny,n1,n2
  330 format(' #no. of extended rows or columns exceeds 2048:'/
     1' input no. of rows and columns='2i4/
     2' no. of rows and columns required for filtering=',2i4,/)
  340 close(11)
      go to 10
  350 write(iw,360)
  360 format(' end of job')
      stop
      end
c***********************************************************************
c    subroutine 'snamel' reads parms from namelist section
c    iread=5 (terminal read) or 10 (read from file10)
c***********************************************************************
      subroutine snamel(iread)
      character*50 fname,coname
      character*2 drive
      common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
      common/parm1/fname,coname,iopt1,iopt2,idval,icoef,ddx,ddy
      common/parm2/den,dec,xinc,bmag,bdec,binc,z,thet1,thet2,istr,w1,w2,
     1w3,w4,drive
c      common/parms/iopt1,iopt2,den,dec,xinc,bmag,bdec,binc,idval,ddx,
c     1ddy,xo,yo,thet1,thet2,istr,nadd,w1,w2,w3,w4,icoef,z
c                               do something about end = 10
        call namemc(iread)
c      read(iread,parms,end=10)
      return
   10 write(iw,20)
   20 format('#odd eof in parms command file--program aborted')
      if(iread.eq.10)close(10)
      stop
      end
























c***********************************************************************
      subroutine spspol(itype,x,y,xysq,ar,ai)
c       subroutine 'spspol' operates on input fourier coefficients to
c       obtain the magnetic field reduced to the pole, to make a
c       pseudo-gravity (or pseudo-magnetic) transformation, or convert
c       vertical-component data to total-field data.
c***********************************************************************
      common/basic/c1,el,em,en,bl,bm,bn
c itype   operator
c   1      psdgrv
c   2      redpol
c   8      psdmag
c  10      vertot
c
      xy=sqrt(xysq)
      if(itype-9) 5,50,45
    5 bml=x*bl+y*bm
      eml=x*el+y*em
      ber=bn*en*xy-bml*eml/xy
      bei=(en*bml+bn*eml)
      go to (10,20,50,50,50,50,50,40), itype
   10 c=c1
      go to 30
   20 c=xy
   30 art=ar
      ait=ai
      bes=c/(ber*ber+bei*bei)
      ar=(ber*art+bei*ait)*bes
      ai=(ber*ait-bei*art)*bes
      go to 50
   40 c=1./c1
      art=ar
      ait=ai
      ar=(ber*art-bei*ait)*c
      ai=(bei*art+ber*ait)*c
      go to 50
c new vertot option, code taken from Hildenbrand's vertot.for
   45 ber=en
      bei=(el*x+em*y)/xy
      art=ar
      ait=ai
      ar=(ber*art-bei*ait)
      ai=(bei*art+ber*ait)
c
   50 return
      end
c***********************************************************************
       subroutine sfvert(c4,xysq,ar,ai)
c  subroutine 'fvert' operates on input fourier coefficients
c  to obtain the 1st vertical derivative of the input data.
c***********************************************************************
       xsqrt=sqrt(xysq)
       xyp=xsqrt*c4
       ar=ar*xyp
       ai=ai*xyp
       return
       end
c***********************************************************************
       subroutine ssvert(c4,xysq,ar,ai)
c  subroutine 'svert' operates on input fourier coefficients
c  to obtain the 2nd vertical derivative of the input data.
c***********************************************************************
       xyp=xysq*c4
       ar=ar*xyp
       ai=ai*xyp
       return
       end
c***********************************************************************
      subroutine sbpass(a1,a2,f1,f2,f3,f4,btan12,btan34,xysq)
c   subroutine 'sbpass' eliminates wavelengths
c   lying outside the range of w1,w2,w3 & w4.
c***********************************************************************
      rad=sqrt(xysq)
      if(rad.lt.f1) go to 10
      if(rad.gt.f4) go to 10
      go to 20
   10 a1=0.0
      a2=0.0
      go to 40
   20 if(rad.ge.f2.and.rad.le.f3) go to 40
      if(rad.gt.f2) go to 30
      f=(rad-f1)*btan12
      a1=f*a1
      a2=f*a2
      go to 40
   30 f=(f4-rad)*btan34
      a1=f*a1
      a2=f*a2
   40 return
      end
c***********************************************************************
      subroutine strend(x,y,ar,ai)
c  subroutine 'trend' passes or rejects trends striking
c  between angles thet1 and thet2.
c***********************************************************************
      common/tr/prad,p,q,th1p,th2p
      if(x.eq.0.0) go to 50
      if(y.eq.0.0) go to 60
      thet=atan(x/y)*prad
      if(thet.lt.0.0)thet=180.+thet
   10 if(th2p.gt.th1p) go to 40
      if(thet.gt.th2p.and.thet.lt.th1p) go to 30
   20 ar=ar*q
      ai=ai*q
      return
   30 ar=ar*p
      ai=ai*p
      return
   40 if(thet.gt.th1p.and.thet.lt.th2p) go to 30
      go to 20
   50 thet=0.
      if(y.lt.0.0)thet=180.
      go to 10
   60 thet=90.
      go to 10
      end
c***********************************************************************
      subroutine srmean(nx,np,a1)
c  subroutine 'srmean' removes the mean from the input grid
c  employing only boundary values.
c***********************************************************************
      dimension a1(np,3),title(14)
      character*50 fname,coname,fileo
      real mean
      common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
      common/parm1/fname,coname,iopt1,iopt2,idval,icoef,ddx,ddy
      data dval/1.0e37/
      ko=10
      mean=0.0
      k=0
      if(idval.eq.0) go to 150
c  idval>0: find and sum boundary values.
      call srdbin(kr,a1,ny)
      do 10 j=1,ny
      if(a1(j,1).ge.dval) go to 10
      mean=mean+a1(j,1)
      k=k+1
   10 continue
      call srdbin(kr,a1(1,2),ny)
      do 130 i=3,nx
      call srdbin(kr,a1(1,3),ny)
      iq=-1
      if(a1(1,2).ge.dval) go to 20
      k=k+1
      mean=mean+a1(1,2)
      if(a1(1,1).lt.dval.and.a1(1,3).lt.dval)iq=0
   20 do 90 j=2,ny-1
      if(a1(j,2).ge.dval) go to 60
      if(iq)30,40,50
   30 mean=mean+a1(j,2)
      k=k+1
      iq=0
      if(a1(j,1).ge.dval.or.a1(j,3).ge.dval)iq=-1
      go to 90
   40 iq=1
      if(a1(j,1).lt.dval.and.a1(j,3).lt.dval) go to 90
      iq=-1
      mean=mean+a1(j,2)
      k=k+1
      go to 90
   50 iq=1
      if(a1(j,1).lt.dval.and.a1(j,3).lt.dval) go to 90
      mean=mean+a1(j,2)+a1(j-1,2)
      k=k+1
      iq=-1
      go to 90
   60 if(iq)80,80,70
   70 mean=mean+a1(j-1,2)
      k=k+1
   80 iq=-1
   90 continue
      if(a1(ny,2).ge.dval) go to 100
      k=k+1
      mean=mean+a1(ny,2)
      if(a1(ny,1).lt.dval.and.a1(ny,3).lt.dval) go to 110
  100 if(iq.ne.1) go to 110
      mean=mean+a1(ny-1,2)
      k=k+1
  110 do 120 j=1,ny
      a1(j,1)=a1(j,2)
  120 a1(j,2)=a1(j,3)
  130 continue
      do 140 j=1,ny
      if(a1(j,2).ge.dval) go to 140
      mean=mean+a1(j,2)
      k=k+1
  140 continue
      go to 240
c  idval=0: read grid and check boundary values for flagged values.
  150 i=0
  160 i=i+1
      call srdbin(kr,a1,ny)
      if(i.gt.1.and.i.lt.nx) go to 190
      do 170 j=1,ny
      if(a1(j,1).ge.dval) go to 180
  170 mean=mean+a1(j,1)
      if(i-nx)160,230,160
  180 write(iw,220)i,j,a1(j,1)
      close(11)
      stop
  190 if(a1(1,1).ge.dval.or.a1(ny,1).ge.dval) go to 200
      mean=mean+a1(1,1)+a1(ny,1)
      go to 160
  200 write(iw,210)i,a1(1,1),a1(ny,1)
  210 format(/,' #found flagged value: row = ',i3,2x,'col = 1 or ny',
     12x,' value = ',2e10.4,' may be more?')
  220 format(/,' #found flagged value: row = ',i3,2x,' col = ',i3,
     12x,'value = ',e10.4' may be more?')
      close(11)
      stop
  230 k=(nx+ny-2)*2
c  compute the mean.
  240 c=mean/float(k)
      ip=iopt1
      if(ip.eq.0)ip=iw
      write(ip,250)c,k
  250 format(/,' mean value subtracted = ',e16.8,5x,i5,' values used',/)
      close(11)
c  read input array and write array with the mean subtacted.
      open(11,status='old',form='unformatted',file=fname,share='denywr')
      read(kr)title,pgm,ny,nx,nz,yop,dyp,xop,dxp
      fileo='fftfil.grd'
      open(10,status='unknown',form='unformatted',file=fileo)
      write(ko)title,pgm,ny,nx,nz,yo,dy,xo,dx
      do 270 i=1,nx
      call srdbin(kr,a1,ny)
      do 260 j=1,ny
  260 a1(j,1)=a1(j,1)-c
  270 call swrbin(ko,a1,ny)
      close(11)
      close(10)
      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)
      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)
c
c       assigns values to proper variables
c       variables are passed to program acdepth 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 fftfil
c
        parameter (nnvar=24,numr=7)
        character*6 pvar,var(nnvar)
        character*56 tvar, kvar, cfmt
        character*50 fname,coname
        character*2 drive
        logical chv  
        common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nadd,ir
        common/parm1/fname,coname,iopt1,iopt2,idval,icoef,ddx,ddy        
        common/parm2/den,dec,xinc,bmag,bdec,binc,z,thet1,thet2,istr,
     1  w1,w2,w3,w4,drive
        data var/'nadd','iopt1','iopt2','idval','icoef','istr','yo',
     1  'xo','ddx','ddy','den','dec','xinc','bmag','bdec','binc',
     2  'z','thet1','thet2','w1','w2','w3','w4','drive'/
        numa= 25
        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 variable kvar
c
        m = 57 - nn
        im = m - 1
        kvar(m:56) = tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write(cfmt, 50) im, nn
50      format ('(',i2,'x,i',i2,')')
        read(kvar, cfmt) jvar
        else
c
c       real value
c
        write (cfmt, 60) im, nn
60      format ('(',i2,'x,g',i2, '.0)')
        read (kvar, cfmt) xvar
        endif
        endif
        goto (101,102,103,1031,104,105,106,107,108,109,110,111,112,113,
     1  114,115,116,117,118,119,120,121,122,123),i
101     nadd = jvar
        goto 200
102     iopt1 = jvar
        goto 200
103     iopt2 = jvar
        goto 200
1031    idval = jvar
        goto 200
104     icoef = jvar
        goto 200
105     istr = jvar
        goto 200
106     yo = xvar
        goto 200
107     xo = xvar
        goto 200
108     ddx = xvar
        goto 200
109     ddy = xvar
        goto 200
110     den = xvar
        goto 200
111     dec = xvar
        goto 200
112     xinc = xvar
        goto 200
113     bmag = xvar
        goto 200
114     bdec = xvar
        goto 200
115     binc = xvar
        goto 200
116     z = xvar
        goto 200
117     thet1 = xvar
        goto 200
118     thet2  = xvar
        goto 200
119     w1 = xvar
        goto 200
120     w2 = xvar
        goto 200 
121     w3 = xvar
        goto 200
122     w4 = xvar
        goto 200
123     drive=tvar(1:1)//':'
        go to 200

190     continue
        write(*,*)' error in namelist - ',pvar,' variables not included'
        stop
200     nvar = i
        return
        end
