c *******************************************************************
c  taylor.for
c
c  Driver for level-to-drape and drape-to-level continuation using 
c  Taylor series expansion.
c
c  V.J.S. Grauch, USGS, March 1984.
c
c     f     =     input data
c     g     =     1st vertical derivative of f
c     h     =     2nd vertical derivative of f
c     z     =     surface data
c     iopt  =     1 if level-to-drape chosen, 2 if drape-to-level
c
c       Be sure to link with increase stack frame size,
c       link /ST:4096 taylor.obj+sfftmg.obj+sfourt.obj+namemc.obj
c *******************************************************************
      character*50 ifile,ofile,outfile,ifilea,ifileb,sfile
      character*56 title
      dimension work(2048),f(2,1024,16),g(2,1024,16),h(2,1024,16),
     1 z(1024),id(14),pgm(2)
      common/l2dparms/flevel,convf,zlevel,ifilt,w1,w2,idval,
     1 nadd,nterms,ifile,ofile,sfile,title,iter
      common/l2dinfo/p1,p2,nz,yo,xo,dx,dy,ncol,nrow
      common/units/iw,ir,in,iout,icmd,is,iunita,iunitb,iunitc,
     1 if,iunitg,iunith
c      data iw/6/,ir/5/,in/10/,iout/12/,icmd/9/,is/11/,iunita/20/,
c     1iunitb/21/,iunitc/22/,if/14/,iunitg/23/,iunith/24/,lnri/16/
        iw=6
        ir=5
        in=10
        iout=12
        icmd=9
        is=11
        iunita=20
        iunitb=21
        iunitc=22
        iunith=24
        if=14
        iunitg=23
        data lnri/16/
1     write(iw,800)
800   format(1x,'1 - level to drape continuation'/' 2 - drape to level
     1 continuation'/8x,'Enter option: ',$)
      read(ir,*) iopt
      if(iopt.lt.1.or.iopt.gt.2) go to 1
c
c  set up input parameters
c
      call setup(iopt)
c  find no. of rows for fft and extend cols by nadd
      call extend(nadd,n1,n2,id2,lnri,nri,nxa)
      if(iopt.eq.2.or.iopt.eq.4) go to 20
c
c   LEVEL TO DRAPE OPERATION
c
      p1='*l t'
      p2='o d*'
      open(iout,file=ofile,form='unformatted',status='unknown')
      write(iout) title,p1,p2,ncol,nrow,nz,xo,dx,yo,dy
      write(iw,804)
804   format(/' Computing...'/)
      call l2d(f,g,h,work,id2,nri,n1,n2,z,nxa)
      close(in)
      close(is)
      close(iout)
      go to 300
c
c  DRAPE TO LEVEL OPERATION
c
20    finlevel=flevel
      flevel=zlevel
      isave=0
c get 1st approx by using level-to-drape with neg. z
      convf=-convf
      if(iter.eq.-1) then
         isave=1
         iter=2
      endif
      open(iout,file='f1.tmp',status='unknown',form='unformatted')
      write(iout) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      write(iw,805)
805   format(/28x,'RECOVERY ERROR MAGNITUDE'/' Iteration',10x,'Max',
     1 15x,'Mean',12x,'Std dev.'/)
      call l2d(f,g,h,work,id2,nri,n1,n2,z,nxa)
      close(in)
      rewind is
      close(iout)
c
c  Get next approximations
c
      convf=-convf
      outfile='f1d.tmp'
      ifilea='f1.tmp'
      k=0
25    k=k+1
      open(in,file=ifilea,status='old',form='unformatted',mode='read')
      read(in) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      read(is) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      open(iout,file=outfile,status='new',form='unformatted')
      write(iout) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      call l2d(f,g,h,work,id2,nri,n1,n2,z,nxa)
      close(in)
      close(iout)
c  find error and add error to one approx to get next approx
      if(k.eq.2) go to 32
      outfile='f2.tmp'
      ifilea='f1d.tmp'
      ifileb='f1.tmp'
      go to 33
32    ifilea='f2d.tmp'
      ifileb='f2.tmp'
      outfile='f3.tmp'
33    if(k.ne.iter) go to 35
      close(is)
      if(abs(finlevel-zlevel).gt.0.0001) go to 36
      open(iout,file=ofile,form='unformatted',status='unknown')
      p1='*d t'
      p2='o l*'
      write(iout) title,p1,p2,ncol,nrow,nz,xo,dx,yo,dy
      go to 40
35    rewind(is)
36    open(iout,file=outfile,form='unformatted',status='unknown')
      write(iout) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
40    open(in,file=ifile,status='old',form='unformatted',mode='read')
      open(iunita,file=ifilea,status='old',form='unformatted')
      open(iunitb,file=ifileb,status='old',form='unformatted')
      read(in) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      read(iunita) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      read(iunitb) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      errmax=-1.e+38
      sumerr=0.0
      sumerr2=0.0
      en=0.0
      do 50 j=1,nrow
      read(in) dum,(f(1,i,1),i=1,ncol)
      read(iunita) dum,(g(1,i,1),i=1,ncol)
      read(iunitb) dum,(h(1,i,1),i=1,ncol)
      do 45 i=1,ncol
      if(f(1,i,1).gt.1.e+38) go to 45
      err=f(1,i,1)-g(1,i,1)
      if(abs(err).gt.errmax) errmax=abs(err)
      sumerr=sumerr+abs(err)
      if(sumerr.ge.1.e+20) go to 299
      sumerr2=sumerr2+err*err
      en=en+1.0
      f(1,i,1)=err+h(1,i,1)
45    continue
      write(iout) dum,(f(1,i,1),i=1,ncol)
50    continue
      sumerr=sumerr/en
      stddev=(sumerr2-(sumerr*sumerr)/en)/(en-1.0)
      stddev=sqrt(stddev)
      write(iw,802) k,errmax,sumerr,stddev
802   format(5x,i1,9x,g14.6,2(5x,g14.6))
      close(in)
      close(iunita,status='delete')
      if(isave.ne.0) close(iunitb)
      if(isave.eq.0) close(iunitb,status='delete')
      close(iout)
      if(k.eq.iter) go to 100
      ifilea=outfile
      outfile='f2d.tmp'
      go to 25
100   if(isave.ne.0) write(iw,803)
803   format(/' 1st-iter. approx. (on z ref. level) is in file f1.tmp')
      if(abs(finlevel-zlevel).le.0.0001) go to 300
c  level-to-level continuation if needed for final level file
      flevel=zlevel
      zlevel=finlevel
      nterms=1
      open(in,file=outfile,status='old',form='unformatted')
      read(in) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      open(iout,file=ofile,status='new',form='unformatted')
      p1='*d t'
      p2='o l*'
      write(iout) title,p1,p2,ncol,nrow,nz,xo,dx,yo,dy
      call l2d(f,g,h,work,id2,nri,n1,n2,z,nxa)
      close(in,status='delete')
      close(iout)
      go to 300
299   write(iw,806)
806   format(' Method is rapidly diverging.  Abort.')
      close(iunita,status='delete')
      close(iunitb,status='delete')
      close(iout,status='delete')
300   stop
      end
c **********************************************************************
c
      subroutine setup(iopt)
c sets up the input parameters for the level to drape or
c drape to level operations.  Used with taylor.for
c
c ********************************************************************
      character answer*1,title*56
      character*50 ifile,ofile,sfile,command
      common/l2dparms/flevel,convf,zlevel,ifilt,w1,w2,
     1idval,nadd,nterms,ifile,ofile,sfile,title,iter
      common/l2dinfo/p1,p2,nz,yo,xo,dx,dy,ncol,nrow
      common/units/iw,ir,in,iout,icmd,is,iunita,iunitb,iunitc,
     1 if,iunitg,iunith
      dimension id(14),pgm(2),z(1024)
c
c
      ifile=' '
      ofile=' '
      sfile=' '
      title='this is the default title'
1     write(iw,800)
800   format(' Command or command file name: ',$)
      read(ir,801) command
801   format(a50)
      if(command.eq.'exit'.or.command.eq.'EXIT') go to 300
      if(command.eq.'help'.or.command.eq.'HELP'.or.command.eq.' ') then
         write(iw,888)
         go to 1
  888 format(' Commands are'/,5x,'term - enter the parameters desired 
     1interactively and create .cmd file at'/,10x,' the same time.',/,
     2 5x,'help'/ 5x,'exit'/' Anything else entered will be taken to
     3 be a command file name.')
      endif
      if(command.ne.'term'.and.command.ne.'TERM') go to 50
c
c  interactively set up parameters and command file.  command file will
c  be named taylor.cmd
      command='taylor.cmd'
      write(iw,802)
802   format(/,' Begin entering parms interactively.',/,' Want me to
     1 ask verbose (as opposed to brief) questions? '$)
      read(ir,803) answer
803   format(a1)
      if(answer.eq.'y'.or.answer.eq.'Y') iopt=iopt+2
c ask for parms
      if(iopt.eq.2.or.iopt.eq.4) write(iw,855)
855   format(/' Will either input file have dvals? '$)
      if(iopt.eq.1.or.iopt.eq.3) write(iw,856)
856   format(/' Will input data file have dvals? ',$)
      read(ir,803) answer
      idval=0
      if(answer.eq.'y'.or.answer.eq.'Y') idval=1
      if(iopt.gt.2) write(iw,709)
709   format(/'  What is the factor to convert the units of measure of
     1 the surface into the',/,'  units of measure of the grid interval
     2 (e.g. .001 if surface is in m and'/'  grid units in km).  Enter
     3 the negative of this number if the sense of the'/'  surface file
     4 is z positive down.')
      write(iw,809)
809   format(' conversion factor = ',$)
      read(ir,*) convf
      if(iopt.eq.2.or.iopt.eq.4) go to 10
      if(iopt.eq.3) write(iw,706)
706   format(/'  Give the level at which the input data were observed in
     1 the SAME units as those'/'  of the draping surface.  (Remember to
     2 be consistent with whether up or down'/'  is the positive
     3 direction.)')
      write(iw,806)
806   format(' flevel= ',$)
      read(ir,*) flevel
10    if(iopt.gt.2) write(iw,711)
711   format(/'  The Taylor series expands about a z reference level
     1 that should best be near'/'  the mean of the draping surface so
     2 that the deviation of the surface'/'  above and below that level
     3 is the least amount.  Enter the level in the'/'  SAME units used
     4 for the draping surface.')
      write(iw,811)
811   format(' Enter zlevel (1.e+38 to use default):  ',$)
      read(ir,*) zlevel
      kflag=0
      if(zlevel.ge.1.e+38) kflag=1
15      if(iopt.eq.1.or.iopt.eq.3) go to 17
      if(iopt.gt.2) write(iw,712)
712   format(/'  Enter the level to which to continue the data in the
     1 SAME units'/'  as the draping surface units.')
      write(iw,812)
812   format(' flevel = ',$)
      read(ir,*) flevel
17      nterms=3
      if(iopt.eq.1.or.iopt.eq.3) go to 18
      if(iopt.gt.2) write(iw,730)
730   format(/'  Read documentation for examples on when to use
     1 2 or 3 terms of the Taylor'/'  series expansion.')
      write(iw,830)
830   format(' Want 2 or 3 terms in Taylor expansion? ',$)
      read(ir,*) nterms
18    if(iopt.gt.2) write(iw,713)
713   format(/'  Bandpass (lowpass) filtering is often desirable during
     1 the 1st & 2nd vert.'/'  derivative-taking of the Taylor series
     2 expansion or in the continuation'/'  of flevel to zlevel or vice
     3 versa.  Of course if zlevel=flevel, filtering'/'  the continuation
     4 is not desirable.')
20    write(iw,813)
813   format(' Enter filtering option (-1 for help): ',$)
      read(ir,*) ifilt
      if(ifilt.eq.0) go to 30
      if(ifilt.gt.0.and.ifilt.le.3) go to 25
      write(iw,833)
833   format('   0 - no filtering'/'   1 - filter derivs only'/'   2 -
     1 filter derivs and continuation'/'   3 - filter continuation
     2 only')
      go to 20
25    if(iopt.gt.2) write(iw,714)
714   format(/'  Two wavelengths w1<=w2 describe a (tapered) lowpass
     1 filter.  Remember',/,'  wavelength=2 * anomaly width.  Defaults:
     2  w1=2*grid interval, w2=w1+2*grid'/'  interval.  w1 &
     3 w2 are in SAME units as grid interval.')
      write(iw,814)
814   format(' Enter w1 and w2 (type 0 0 to use default) ',$)
      read(ir,*) w1,w2
      if(w1.gt.w2) go to 25
      if(w2.ne.0.) go to 30
      w1=1.e+38
      w2=1.e+38
      kflag=kflag+2
30    if(iopt.gt.2) write(iw,724)
724   format(/'  It often helps the FFT operation to extend the number
     2 of rows and',/,'  columns by a small proportion of the total grid.
     3  Enter the number to'/'  add to the cols and rows:')
      write(iw,824)
824   format(' nadd = ',$)
      read(ir,*) nadd
      if(iopt.eq.1.or.iopt.eq.3) go to 32
31    if(iopt.gt.2) write(iw,725)
725   format(/'  The approximation of data on the level is found through
     1 successive'/'  iterations (1 or 2).  More iterations doesn''t
     2 neccessarily mean a better'/'  approximation though (see docu
     3mentation).  Enter iter=-1 if'/'  you want 2 iterations but also
     4 wish to save the grid of the 1st iteration.')
      write(iw,825)
825   format(' Enter iter: ',$)
      read(ir,*) iter
      if(iter.eq.-1.or.iter.le.2.and.iter.ge.1) go to 32
      write(iw,725)
      go to 31
c
c  create command file if requested
c
32    if(iopt.gt.2) write(iw,715)
715   format(/' Now I can create a command file with all the
     1 parameters you picked',/,' so that you can edit and/or use it
     2 for a later program run.')
      write(iw,860)
860   format(/' Want to save this info in command file? ',$)
      read(ir,803) answer
      if(answer.ne.'y'.and.answer.ne.'Y') go to 51
      open(icmd,file=command,status='unknown',form='formatted')             
      if(kflag.gt.0) go to 35
      write(icmd,818) flevel,convf,nterms,zlevel,ifilt,w1,w2,nadd,idval
818   format(' &parms  flevel=',g14.6,',convf=',g14.6,',nterms=',i1,
     2',zlevel=',g14.6,',',/,' ifilt=',i1,',w1=',g14.6,',w2=',g14.6,/,
     3' nadd=',i2,',idval=',i1,' &')
      go to 40
35    if(kflag.gt.1) go to 37
      write(icmd,819) flevel,convf,nterms,ifilt,w1,w2,nadd,idval
819   format(' &parms flevel=',g14.6,',convf=',g14.6,',nterms=',i1,
     1/,' ifilt=',i1,',w1=',g14.6,',w2=',g14.6,/,' nadd=',i2,',idval=',
     2 i1,' &')
      go to 40
37    if(kflag.eq.2) go to 38
      write(icmd,719) flevel,convf,nterms,ifilt,nadd,idval
719   format(' &parms flevel=',g14.6,',convf=',g14.6,',nterms=',i1,
     1/,' ifilt=',i1,',nadd=',i2,',idval=',i1,' &')
      go to 40
38    write(icmd,720) flevel,convf,nterms,zlevel,ifilt,nadd,idval
720   format(' &parms flevel=',g14.6,',convf=',g14.6,',nterms=',i1,
     1/' zlevel=',g14.6,',ifilt=',i1,',nadd=',i2,',idval=',i1,' &')
40    close(icmd)
      write(iw,820)
820   format(/' **command file taylor.cmd created**')
      go to 51
c
c  get parms from existing command file
c
c  set default values of parms
50    nadd=0
      convf=1.
      ifilt=0
      flevel=0.
      zlevel=1.e+38
      w1=1.e+38
      w2=1.e+38
      nterms=3
      iter=-3
      open(icmd,file=command,form='formatted',status='old',
     1mode='read',err=299)
      call namemc(icmd)
c check for parameter errors
      if(ifilt.lt.0.or.ifilt.gt.3) go to 290
      if(ifilt.ne.0.and.w1.gt.w2) go to 290
      if(idval.ne.0.and.idval.ne.1) go to 290
      if(nadd.lt.0) go to 290
      if(nterms.lt.1.or.nterms.gt.3) go to 290
      if(iopt.eq.1.or.iopt.eq.3) go to 51
      if(iter.eq.-1.or.iter.eq.2.or.iter.eq.1) go to 51
      write(iw,825)
      read(ir,*) iter
49    if(iter.eq.-1.or.iter.eq.2.or.iter.eq.1) go to 51
      write(iw,725)
      write(iw,825)
      read(ir,*) iter
      go to 49
c prompt for input files
51    if(ifile.ne.' ') go to 55
52    if(iopt.eq.1.or.iopt.eq.3) write(iw,704)
704   format(/' Input file (level survey): ',$)
      if(iopt.eq.2.or.iopt.eq.4) write(iw,804)
804   format(/' Input file (draped survey): ',$)
      read(ir,801) ifile
55    open(in,file=ifile,status='old',form='unformatted',mode='read',
     1err=297)
      read(in)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      if(dx.ne.dy) go to 295
      if(sfile.ne.' ') go to 57
56    write(iw,807)
807   format(/' Enter draping or irregular surface grid: ',$)
      read(ir,801) sfile
57    open(is,file=sfile,form='unformatted',mode='read',status='old',
     1err=298)
      read(is) id,pgm,ncol2,nrow2,nz,xo2,dx2,yo2,dy2
      if(abs(xo-xo2).gt.0.001.or.abs(yo-yo2).gt.0.001) go to 293
      if(abs(dx-dx2).gt.0.001*abs(dx).or.dx2.ne.dy2) go to 293
      if(nrow.ne.nrow2.or.ncol.ne.ncol2) go to 293
      if(ofile.ne.' ') go to 59
      write(iw,815)
815   format(/' Output file: ',$)
      read(ir,801)ofile
59    if(title.ne.'this is the default title') go to 100
      write(iw,816)
816   format(' Enter title')
      read(ir,817) title
817   format(a56)
c check zero level or set it if default used  
100   call reflvl(is,reclevel,nrow,ncol,dx,iw,convf,z)
      if(zlevel.eq.1.e+38) go to 105
      if(abs((zlevel-reclevel)*convf).gt.abs(1.5*dx)) write(iw,821)
     1 zlevel,reclevel
821   format(/' **WARNING: zero level of ',g14.6,' is more than 1.5*dx
     1 away from',/,'    recommended level of ',g14.6)
      go to 110
105   zlevel=reclevel
      write(iw,822) zlevel
822   format(/' Using z ref. level of ',g14.6)
c  set w1 and w2 if default values picked
110   if(w1.lt.1.e+38.and.w2.lt.1.e+38.or.ifilt.eq.0) go to 300
      w1=2*dx
      w2=4*dx
      write(iw,823) w1,w2
823   format(' Using w1,w2= ',g14.6,1x,g14.6)
      go to 300
c
c ERROR  MESSAGES
c
290   write(iw,292)
292   format(' Command file has parameter error.  Abort.')
      stop
293   write(iw,294)
294   format(' Grids don''t match')
      stop
295   write(iw,296)
296   format(' dx must = dy')
      stop
297   write(iw,291) ifile
291   format(3x,a12,' not found or not binary.  Try again. ')
      go to 52
298   write(iw,291) sfile
      go to 56
299   write(iw,391) command
391   format(' File not found or not ascii.  Try again.')
      go to 1
300   return
      end
c ********************************************************************
c      This subroutine recommends a reference level for the 
c       level-to-drape operation.
c **********************************************************************
      subroutine reflvl(is,rlevel,nrow,ncol,dx,iw,convf,z)
      dimension z(ncol),id(14),pgm(2)
      data zmax/-1.e38/,zmin/1.e38/
c
c  calculate min, max and mean of the draping surface
c
      sumz=0.0
      en=0.0
      do 20 j=1,nrow
      read(11,end=80,err=70) dum,z
      do 10 i=1,ncol
c convert surf. file to z positive down and into units of dx
c calculate mean and std dev
      if(z(i).ge.1.e+38) go to 10
      z(i)=z(i)*convf
      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
      if(j.eq.1.and.en.eq.0.0) go to 60
20    continue
      rewind is
      read(is) id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      zmean=sumz/en
c
c  Ideally zmax-rlevel should be <= .9*dx (distance method will need to
c  approximate downward) and rlevel-zmin <= 1.5*dx (distance upward).
c  Try to get close to this rlevel.  Use the mean if conditions are
c  satisfied with rlevel=zmean.  Round off rlevel.
c
      rlevel=zmean
      if(zmax-rlevel.le.0.9*dx.and.rlevel-zmin.le.1.5*dx) go to 50
      rtest1=zmax-.9*dx
      rtest2=zmin+1.5*dx
      rlevel=(rtest1+rtest2)*.5
      if(zmax-rlevel.le.1.2*dx.and.rlevel-zmin.le.2.0*dx) go to 50
      write(iw,800)
800   format(/' **WARNING: Relief of draping surface is greater than
     1 dx.',/,' This could cause the method to diverge.  Try
     2 a coarser dx?')
c  round off rlevel
50    if(rlevel.eq.0.0) go to 55
      nsign=rlevel/abs(rlevel)
      ntest=2-ifix(alog10(abs(rlevel)))
      test=10.e+0**float(ntest)
      rtest=float(nint(rlevel*test))
      rlevel=nsign*rtest/test
c change units of rlevel back and return
55    rlevel=rlevel/convf
      return
60    write(iw,801)
801   format(' No data in entire first row of surface file.  Can''t
     1 handle it.')
      stop
70    write(iw,802)
802   format(' Error in read of surface file.')
      stop
80    write(is,803)
803   format(' unexpected eof in surface file.')
      stop
      end
C***********************************************************************
C  subroutine 'l2d' removes flagged values from grid, transforms
C  data, calls routines to apply filter,  performs level-to-drape
C  Taylor's series expansion, and writes resulting grid.
C  Note that both input data files and output file must be open and
C  have had headers read before calling this subroutine.  They will
C  remain open upon return.
C  Guts of this subroutine lifted from program FFTFIL by Thomas G.
C  Hildenbrand, 1983, USGS Open-File Report 83-237.
C     f     =     input data
C     g     =     storage for 1st vertical derivative of f
C     h     =     storage for 2nd vertical derivative of f
C     work  =     work array for subroutine sfftmg
C     n1    =     number of rows
C     n2    =     number of cols
C     id2   =     2 * max(n1,n2)
C     nri   =
C     z     =     surface data
C     nxa   =
C
C***********************************************************************
      Subroutine l2d(f[huge],g[huge],h[huge],work,id2,nri,n1,n2,z,nxa)
      Dimension f(2,n2,nri),g(2,n2,nri),h(2,n2,nri),work(id2),
     1jflag(2,20),z(n2)
      Character*50 ifile,ofile,sfile
      Character*56 title
      Common /l2dparms/flevel,convf,rlevel,ifilt,w1,w2,idval,nadd,
     1nterms,ifile,ofile,sfile,title,iter
      Common /l2dinfo/p1,p2,nz,yo,xo,dx,dy,ny,nx
      Common /units/iw,ir,in,iout,icmd,is,iunita,iunitb,iunitc,if,
     1iunitg,iunith
      Data pi2/6.28318531/,dval/1.701412E38/
      n22=2*n2
      If (ifilt.EQ.0) Go To 20
      If (w1.NE.0.0) Go To 5
      f4=1.0e+30
      Go To 7
    5 f4=1./w1
    7 If (w2.NE.0.0) Go To 8
      f3=1.0e+30
      Go To 10
    8 f3=1./w2
   10 If (f3.NE.f4) denom=1./(f4-f3)
C  compute beginning and end column & row in original grid.
   20 Do 30 i=1,2
      Do 30 j=1,n2
   30 f(2,j,i)=0.0
      nap1=nadd+1
      n2a=n2-nadd
      nxap1=nxa*0.5+1
      fxap1=float(nxa)/2.+1.
      If (abs(fxap1-float(nxap1)).GT.0.0001)nxap1=nxap1+1
      n1a=nxap1-1+nx
C  remove flagged values and add '2*nadd' rows & columns to
C  reduce the effect of gibbs phenomena.
      Open (iunita,access='direct',status='unknown',form='unformatted',
     1file='slave1.tmp',recl=n22*4)            
      If (idval.NE.0) Open (if,status='unknown',form='unformatted',
     1file='flag.loc')
      irr=1
      iwr=0
      jr=0
   40 jr=jr+1
      If (idval.EQ.0) Go To 310
C   idval=1: removal of flagged values, create file flag.loc
      If (jr.GT.1) Go To 140
C   read 1st row & remove flagged values.
      Read (in)dum,(f(1,i,1),i=nap1,n2a)
      Read (in)dum,(f(1,i,2),i=nap1,n2a)
      nflag=1
      i=nap1
   50 jflag(1,1)=nap1
      If (f(1,nap1,1).LE.1.0e38)Go To 90
      Do 60 i=nap1,n2a
   60 If (f(1,i,1).LE.1.0e38)Go To 70
      Go To 1270
   70 jflag(1,1)=i
      aval=f(1,i,1)
      Do 80 ii=nap1,i-1
   80 f(1,ii,1)=aval
   90 If (i.GE.n2a)Go To 101
      Do 100 ii=i+1,n2a
      If (f(1,ii,1).LE.1.0e38)Go To 100
      jflag(2,nflag)=ii-1
      Go To 110
  100 Continue 
101     continue
      jflag(2,nflag)=n2a
      Write (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      Go To 350
  110 aval=f(1,ii-1,1)
      Do 120 i=ii,n2a
      If (f(1,i,1).LE.1.0e38)Go To 130
  120 f(1,i,1)=aval
      Write (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      Go To 350
  130 f(1,i-1,1)=(f(1,i-1,1)+f(1,i,1))*0.5
      nflag=nflag+1
      jflag(1,nflag)=i
C   1st row completed: now extend grid & write to disk.
      Go To 90
C   remove flagged values from remaining rows.
  140 If (jr.GE.nx)Go To 150
      Read (in)dum,(f(1,i,3),i=nap1,n2a)
      Go To 170
  150 Do 160 i=nap1,n2a
  160 f(1,i,3)=dval
  170 nflag=1
      i=nap1
      jflag(1,1)=nap1
      If (f(1,nap1,2).LE.1.0e38)Go To 220
      Do 180 i=nap1,n2a
  180 If (f(1,i,2).LE.1.0e38)Go To 190
      i=nap1
      jflag(1,1)=n2a+1
      If (f(1,nap1,3).GE.1.0e38)f(1,nap1,2)=f(1,nap1,1)
      If (f(1,nap1,3).LE.1.0e38)f(1,nap1,2)=(f(1,nap1,1)+
     1f(1,nap1,3))*0.5
      Go To 220
  190 jflag(1,1)=i
      ii=i
      Do 210 irev=nap1,i-1
      ii=ii-1
      If (f(1,ii,3).GE.1.0e38)Go To 200
      f(1,ii,2)=(f(1,ii,1)+f(1,ii+1,2)+f(1,ii,3))*0.3333333
      Go To 210
  200 f(1,ii,2)=(f(1,ii,1)+f(1,ii+1,2))*0.5
  210 Continue
  220 If (i.GE.n2a)Go To 231
      Do 230 ii=i+1,n2a
      If (f(1,ii,2).LE.1.0e38)Go To 230
      jflag(2,nflag)=ii-1
      Go To 240
  230 Continue  
231     continue
      jflag(2,nflag)=n2a
      Write (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      Go To 390
  240 Do 290 i=ii,n2a
      If (i.EQ.n2a)Go To 250
      If (f(1,i+1,2).LE.1.0e38)Go To 270
  250 If (f(1,i,3).GE.1.0e38)Go To 260
      f(1,i,2)=(f(1,i,1)+f(1,i-1,2)+f(1,i,3))*0.3333333
      Go To 290
  260 f(1,i,2)=(f(1,i,1)+f(1,i-1,2))*0.5
      Go To 290
  270 If (f(1,i,3).GE.1.0e38)Go To 280
      f(1,i,2)=(f(1,i,1)+f(1,i-1,2)+f(1,i+1,2)+f(1,i,3))*0.25
      Go To 300
  280 f(1,i,2)=(f(1,i,1)+f(1,i-1,2)+f(1,i+1,2))*0.3333333
      Go To 300
  290 Continue
      Write (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      Go To 390
  300 nflag=nflag+1
      i=i+1
      jflag(1,nflag)=i
      Go To 220
C   idval.eq.0:  check input data.
  310 If (jr.NE.1)irr=2
      Read (in)dum,(f(1,i,irr),i=nap1,n2a)
      Do 320 ii=nap1,n2a
  320 If (f(1,ii,irr).GE.1.0e38)Go To 330
      If (jr.NE.1)Go To 390
      Go To 350
  330 Write (iw,340)jr,ii
  340 Format (' #flagged value row=',i4,3x,'col=',i4,' may be more?')
      Close (in)
      Close (iunita)
      Return
C   extend row & write to disk.
  350 Do 360 i=1,nap1
  360 f(1,i,1)=f(1,nap1,1)
      Do 370 i=n2a,n2
  370 f(1,i,1)=f(1,n2a,1)
  380 iwr=iwr+1
      Call swrda(iunita,iwr,f,n22)
      If (iwr.EQ.nxap1)Go To 40
      Go To 380
  390 Do 400 i=1,nap1
  400 f(1,i,2)=f(1,nap1,2)
      Do 410 i=n2a,n2
  410 f(1,i,2)=f(1,n2a,2)
      iwr=iwr+1
      Call swrda(iunita,iwr,f(1,1,2),n22)
      If (iwr.GE.n1a)Go To 430
      If (idval.LT.1)Go To 40
      Do 420 i=nap1,n2a
      f(1,i,1)=f(1,i,2)
  420 f(1,i,2)=f(1,i,3)
      Go To 40
  430 If (iwr.GE.n1)Go To 450
  440 iwr=iwr+1
      Call swrda(iunita,iwr,f(1,1,2),n22)
      If (iwr.LT.n1)Go To 440
  450 If (idval.NE.0)Close (if)
  460 Format (i5)
C  obtain complex fourier coefficients (f.f.t. f).
  590 Open (iunitb,access='direct',status='unknown',form='unformatted',
     1file='slave2.tmp',recl=n22*4)
C
C    Transform data
C
  980 isign=-1
      Call sfftmg(iunita,iunitb,n2,n1,f,nri,isign,work)
C  f array now holds complex fourier coefficients.
C  start filtering
 1000 Go To (1005,1003,1001),nterms
 1001 Open (iunith,access='direct',status='unknown',form='unformatted',
     1file='slave4.tmp',recl=n22*4)
 1003 Open (iunitg,access='direct',status='unknown',form='binary',
     1file='slave3.tmp',recl=n22*4)
 1005 d=(flevel-rlevel)*convf
      cd=pi2*d
      c1=pi2
      n1p1=n1+1
      n2p1=n2+1
      nnh1=(n1/2)+1
      nnh2=(n2/2)+1
C  operating on f fourier coefficients.
      rn1=1./(float(n1)*dx)
      rn2=1./(float(n2)*dy)
      jr=0
      jmr=n1p1
      igo=1
      If (ifilt.EQ.0)Go To 1010
      igo=2
      If (nterms.EQ.1.AND.ifilt.LE.1)igo=3
      If (nterms.EQ.1.AND.ifilt.GE.2)igo=4
      dumc=1.0
      dumd=1.0
      If (ifilt.LE.2)dumd=0.0
      If (ifilt.GE.2)dumc=0.0
 1010 jr=jr+1
      If (jr.GT.nnh1)Go To 1110
      jj=jr-1
      x=float(jj)*rn1
      xsq=x*x
      Call srdda(iunita,jr,f,n22)                       
      Do 1090 i=1,n2
      ii=i-1
      If (i.GT.nnh2)ii=-(n2p1-i)
      y=float(ii)*rn2
      xysq=xsq+y*y
      w=sqrt(xysq)
      fr=f(1,i,1)
      fi=f(2,i,1)
      Go To (1020,1030,1060,1070),igo
C modify coefs, no bandpass filtering
 1020 e=exp(cd*w)
      e1=e*c1*w
      e2=e1*c1*w
      f(1,i,1)=fr*e
      f(2,i,1)=fi*e
      g(1,i,1)=fr*e1
      g(2,i,1)=fi*e1
      h(1,i,1)=fr*e2
      h(2,i,1)=fi*e2
      Go To 1090
C  modify coefs while bandpass filtering, dumd and dumc are switches to
C  filter or not to filter depending on ifilt
C     dumd=1.0, no filter of derivs; dumc=1.0, no filter of continuation
 1030 If (w.LE.f4)Go To 1050
      e=exp(cd*w)
      f(1,i,1)=fr*e*dumc
      f(2,i,1)=fi*e*dumc
      e1=e*c1*w*dumd
      e2=e1*c1*w
      g(1,i,1)=fr*e1
      g(2,i,1)=fi*e1
      h(1,i,1)=fr*e2
      h(2,i,1)=fi*e2
      Go To 1090
 1050 factor=1.0 
      If (w.GT.f3)factor=(f4-w)*denom
      e=exp(cd*w)
      factrd=(1.0-dumd)*factor+dumd
      factrc=(1.-dumc)*factor+dumc
      e1=e*c1*w*factrd
      e2=e1*c1*w   
      f(1,i,1)=fr*e*factrc
      f(2,i,1)=fi*e*factrc
      g(1,i,1)=fr*e1
      g(2,i,1)=fi*e1
      h(1,i,1)=fr*e2
      h(2,i,1)=fi*e2
      Go To 1090
C continuation only (drape-to-level), no bandpass filtering
 1060 e=exp(cd*w)
      f(1,i,1)=fr*e
      f(2,i,1)=fi*e
      Go To 1090
C  continuation only with bandpass filtering
 1070 If (w.LE.f4)Go To 1075
      e=exp(cd*w)
      f(1,i,1)=0.0
      f(2,i,1)=0.0
      Go To 1090
 1075 factor=1.0
      If (w.GT.f3)factor=(f4-w)*denom
      e=exp(cd*w)*factor
      f(1,i,1)=fr*e
      f(2,i,1)=fi*e
 1090 Continue
      Call swrda(iunita,jr,f,n22)
      Go To (1095,1093,1092),nterms
 1092 Call swrda(iunith,jr,h,n22)
 1093 Call swrda(iunitg,jr,g,n22)
 1095 If (jr.EQ.1)Go To 1010
      jmr=jmr-1
      If (jr.EQ.jmr)Go To 1110
      f(2,1,1)=-f(2,1,1)
      g(2,1,1)=-g(2,1,1)
      h(2,1,1)=-h(2,1,1)
      iwr=n2p1
      Do 1100 i=2,nnh2
      iwr=iwr-1
      x=f(2,i,1)
      y=f(1,i,1)
      f(2,i,1)=-f(2,iwr,1)
      f(1,i,1)=f(1,iwr,1)
      f(1,iwr,1)=y
      f(2,iwr,1)=-x
      x=g(2,i,1)
      y=g(1,i,1)
      g(2,i,1)=-g(2,iwr,1)
      g(1,i,1)=g(1,iwr,1)
      g(1,iwr,1)=y
      g(2,iwr,1)=-x
      x=h(2,i,1)
      y=h(1,i,1)
      h(2,i,1)=-h(2,iwr,1)
      h(1,i,1)=h(1,iwr,1)
      h(1,iwr,1)=y
      h(2,iwr,1)=-x
 1100 Continue
      Call swrda(iunita,jmr,f,n22)
      Go To (1010,1105,1102),nterms
 1102 Call swrda(iunith,jmr,h,n22)
 1105 Call swrda(iunitg,jmr,g,n22)
      Go To 1010
C  f array now holds complex coefficients operated on by filter
C  compute desired anomaly = inv. f.f.t. of f
 1110 isign=1
      Call sfftmg(iunita,iunitb,n2,n1,f,nri,isign,work)
      Go To (1120,1115,1113),nterms            
1113    Call sfftmg(iunith,iunitb,n2,n1,h,nri,isign,work)
1115       Call sfftmg(iunitg,iunitb,n2,n1,g,nri,isign,work)
C  f array now holds computed anomaly on datum = z.
C  do level to drape operation and write to output grid
 1120 If (idval.NE.0)Open (if,status='old',form=
     1'unformatted',file='flag.loc')
      area=1./float(n1*n2)
      dum=0.0
      iwr=0
      dumh=1.0
      If (nterms.EQ.2)dumh=0.0
      If (nterms.EQ.1)Go To 1215
      Do 1210 jr=nxap1,n1a
      iwr=iwr+1
      Call srdda(iunita,jr,f,n22)
      Go To (300,1124,1123),nterms
 1123 Call srdda(iunith,jr,h,n22)
 1124 Call srdda(iunitg,jr,g,n22)
      Read (is)dum,(z(k),k=1,ny)
      If (idval.EQ.0)Go To 1170
      Read (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      j=1
      ik=0
C   restore shape of original grid (i.e. insert flagged values).
      Do 1160 i=nap1,n2a
      ik=ik+1
      If (i.LT.jflag(1,j))Go To 1130
      If (i.EQ.jflag(2,j))Go To 1140
      Go To 1150
 1130 f(1,i,1)=dval
      Go To 1160
 1140 j=j+1
      If (j.LE.nflag)Go To 1150
      j=j-1
      jflag(1,j)=1024
C  convert z(ik) to units of grid interval, positive down, and
C  into units of distance from reference level
 1150 If (z(ik).GE.1.e+38)Go To 1155
      z(ik)=(rlevel-z(ik))*convf
C  Taylor's series expansion.  dumh dummies out the 2nd deriv term
C    depending on nterms
      f(1,i,1)=f(1,i,1)+z(ik)*g(1,i,1)+z(ik)*z(ik)*h(1,i,1)*dumh
      f(1,i,1)=f(1,i,1)*area
      Go To 1160
 1155 f(1,i,1)=dval
 1160 Continue
      Go To 1190
 1170 ik=0
      Do 1180 i=nap1,n2a
      ik=ik+1
      If (z(ik).GE.1.e+38)Go To 1175
      z(ik)=(rlevel-z(ik))*convf
      f(1,i,1)=f(1,i,1)+z(ik)*g(1,i,1)+z(ik)*z(ik)*h(1,i,1)*dumh
      f(1,i,1)=f(1,i,1)*area
      Go To 1180
 1175 f(1,i,1)=dval
 1180 Continue
 1190 Write (iout)dum,(f(1,i,1),i=nap1,n2a)
 1210 Continue
      Go To 1250
C restore shape of grid- continuation only operation
 1215 Do 1247 jr=nxap1,n1a
      iwr=iwr+1
      Call srdda(iunita,jr,f,n22)
      If (idval.EQ.0)Go To 1235
      Read (if)nflag,(jflag(1,j),jflag(2,j),j=1,nflag)
      j=1
      Do 1230 i=nap1,n2a
      If (i.LT.jflag(1,j))Go To 1220
      If (i.EQ.jflag(2,j))Go To 1225
      Go To 1227
 1220 f(1,i,1)=dval
      Go To 1230
 1225 j=j+1
      If (j.LE.nflag)Go To 1227
      j=j-1
      jflag(1,j)=1024
 1227 f(1,i,1)=f(1,i,1)*area
 1230 Continue
      Go To 1245
 1235 Do 1240 i=nap1,n2a
      f(1,i,1)=f(1,i,1)*area
 1240 Continue
 1245 Write (iout)dum,(f(1,i,1),i=nap1,n2a)
 1247 Continue
 1250 If (idval.EQ.1)Close (if,status='delete')
      Close (iunitc,status='delete')
      Close (iunita,status='delete')
      Close (iunitb,status='delete')
      Go To (1290,1265,1263),nterms
 1263 Close (iunith,status='delete')
 1265 Close (iunitg,status='delete')
      Go To 1290
 1270 Write (iw,1280)
 1280 Format (' #entire first row has no data--program aborted')
      Close (in)
      Close (iunita,status='delete')
      If (idval.GT.0)Close (if,status='delete')
 1290 Return
      End
C ********************************************************************
      Subroutine extend(nadd,n1,n2,id2,lnri,nri,nxa)
C  this subroutine extends the cols of grid if nadd not 0 and finds the
C  best no. of rows to use in the fft routines.  Lifted from program
C  FFTFIL by Thomas G. Hildenbrand.
C **********************************************************************
      Common /l2dinfo/p1,p2,nz,yo,xo,dx,dy,ncol,nrow
      n1=nrow
      n2=ncol
      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.1024)Go To 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 n2=n2+2*nadd
      If (n1.GT.1024.OR.n2.GT.1024)Go To 320
      nxa=nxa+2*nadd
      id2=n1
      If (n2.GT.n1)id2=n2
      id2=2*id2
      Go To 300
  320 Write (6,330)nrow,ncol,n1,n2
  330 Format (' no. of extended rows or cols exceeds 1024:',/,' input
     1no. of rows and cols=',2i4,/,'no. of rows and cols required for
     2 filtering=',2i4,/)
  340 Stop
  300 Return
      End
c ********************************************************************
      subroutine srdda(no,ipos,dat[huge],n)
c   this subroutine reads keyed sequential files.
c ********************************************************************
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
c ********************************************************************
      subroutine swrda(no,ipos,dat[huge],n)
c   this subroutine writes keyed sequential files
c ********************************************************************
      dimension dat(n)
      write(no,rec=ipos)dat
      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 variable
c       variables are passed to program Taylor thru common block
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 Taylor
c
        parameter (nnvar=14, numr=6)
        character*6 pvar,var(nnvar)
        character*56 title,tvar,kvar,cfmt
        character*50  ifile,ofile,sfile
        logical chv
        common/l2dparms/flevel,convf,zlevel,ifilt,w1,w2,
     1  idval,nadd,nterms,ifile,ofile,sfile,title,iter
        data var/'ifilt','idval','nadd','nterms','iter','flevel',
     1  'convf','zlevel','w1','w2','ifile','ofile','sfile','title'/
        numa=15
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 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,104,105,106,107,108,109,110,111,112,
     1  113,114),i
101     ifilt = jvar
        goto 200
102     idval = jvar
        goto 200
103     nadd = jvar
        goto 200
104     nterms=jvar
        goto 200
105     iter = jvar
        goto 200
106     flevel = xvar
        goto 200
107     convf = xvar
        goto 200
108     zlevel = xvar
        goto 200
109     w1 = xvar
        goto 200
110     w2 = xvar
        goto 200
111     ifile = tvar(1:nn)
        goto 200
112     ofile = tvar(1:nn)
        goto 200
113     sfile = tvar(1:nn)
        goto 200
114     title = tvar(1:nn)
        goto 200
190     continue
        write (*,*) ' error in namelist - ',pvar,
     1  ' variable not included'
        stop
200     nvar = i
        return
        end


