
      include 'flib.fi'
      include 'fgraph.fi'
      implicit integer*2 (a-z)
      include 'fgraph.fd'
      integer*4 dum1,bgd1
      character*8 buff,passwrd
      data bgd1/1/,fgd1/14/,passwrd/'sB8yr0hF'/
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      call getarg(1,buff,length)
      if(buff.ne.passwrd) then
        print *,'program can only be run from program GEOCON'
        stop
      endif
      call clearscreen($GCLEARSCREEN)
c*****WRITE MAIN FUNCTION KEY MENU
      call funkey
c*****CHECK FOR FUNCTION KEY ACTIVATION
10    call funact
      go to 10
      stop
      end
c***************************************
      subroutine funkey
      implicit integer*2 (a-z)
      parameter(num=8)
      include 'fgraph.fd'
      record/rccoord/curpos
      character* 52 key(num)
      integer*4 dum1,bgd,bgd1,bgd2
      data bgd/15/,bgd2/4/,fgd/0/
      data key/'         PROGRAM  M C G R D',
     & '         minimum curvature gridding program ',
     & ' ',
     & '         F1  = HELP',
     & '         F2  = INTERACTIVE INPUT/EDITING',
     & '         F3  = COMMAND FILE INPUT',
     & '         F4  = RUN GRIDDING PROGRAM',
     & '         ESC = QUIT PROGRAM'/

      bgd1= getbkcolor()
      fgd1=gettextcolor()
      dum1= setbkcolor(bgd)
      dum2=settextcolor(fgd)

      row=1
      do 10 i=1,num-1
      call settextposition(row,4,curpos)
      call outtext(key(i))
      row=row+2
10    continue
      dum1= setbkcolor(bgd2)
      dum2=settextcolor(bgd)
      call settextposition(25,4,curpos)
      call outtext(key(num))

      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)

      return
      end

c***************************************
      subroutine funact
      implicit integer*2 (a-z)
      include 'fgraph.fd'
      integer*4 ierr,ic,ii,bgd1,dum1
      data bgd1/1/,fgd1/14/
      data ic/50/,ii/10/
      data esc/27/

c*****GET SCAN CODE FOR FUNCTION KEYS
      ch=0
      call getchr(ch)
      if(ch.gt.58.and.ch.lt.69) then
      ich = ch-58
      go to (10,20,30,40) ich
      else if(ch.eq.esc) then
      go to 55
      endif
      go to 60

c*****HELP FILE
10    call page('set')
      call helpm
      go to 50

c*****INTERACTIVE INPUT
20    call page('set')
      call edit(ierr)
      if(ierr.ne.0) pause 'press ENTER to return to menu'
      dum = setvideomode($DEFAULTMODE)
      call clearscreen($GCLEARSCREEN)
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      call funkey
      go to 60

c*****USE COMMAND FILE
30    call page('set')
      call input(ic,ic,ierr)
      if(ierr.ne.0) pause 'press ENTER to return to menu'
      go to 50

c*****RUN PROGRAM
40    call page('set')
      call input(ii,ii,ierr)
      if(ierr.ne.0)then
         pause 'press ENTER to return to menu'
         go to 50
      endif
      call mincs(ierr)
      if(ierr.ne.0) pause 'press ENTER to return to menu'

50    call page('res')
      go to 60

c*****QUIT PROGRAM
55    stop

60    return
      end
c***************************************
      subroutine getinp
      include 'flib.fd'
      record/regs$info/regs
      integer*2 ich,icnt
      character*(*) str
      character*1 cr
      data cr/#0d/
      entry getchr(ich)
c*****CHECK FOR KEYBOARD INPUT
c*****RETURNS IF NO KEYBOARD ENTRY
      regs.bregs.ah=11
      call intdosqq(regs,regs)
      if(regs.bregs.al.eq.0) return
      regs.bregs.ah=8
      go to 5
c*****- THIS ENTRY CLEARS BUFFER AND WAITS FOR KEYBOARD ENTRY 
c*****NO ECHO
      entry getch(ich)
      regs.bregs.ah=12
      regs.bregs.al=8
5     call intdosqq(regs,regs)
      if(regs.bregs.al.eq.0) then
c*****EXTENDED KEYBOARD INPUT(E.G., FUNCTION KEYS)
      call intdosqq(regs,regs)
      ich=regs.bregs.al
      else
c*****ASCII CHARACTER INPUT
      ich=regs.bregs.al
      endif
      return
      entry getstr(str,icnt)
c*****CHECK FOR KEYBOARD INPUT
      regs.bregs.ah=11
      call intdosqq(regs,regs)
      if(regs.bregs.al.eq.0) return
c*****GET STRING OF CHARACTERS 
      icnt=0
      do 10 while(str(icnt:icnt).ne.cr)
        regs.bregs.ah=8
        call intdosqq(regs,regs)
        if(regs.bregs.al.eq.0) then
c*****EXTENDED KEYBOARD INPUT(E.G., ARROW KEYS)
          call intdosqq(regs,regs)
          ich=regs.bregs.al
c*****UP ARROW - RETURN
          if(ich.eq.72) then
            icnt = -72
            return
          endif
c*****DOWN ARROW - RETURN
          if(ich.eq.80) then
            icnt = -80
            return
          endif
c*****LEFT ARROW - ECHO A NON-DESTRUCTIVE BACKSPACE
          if(ich.eq.75) then
            icnt=icnt-1
            if(icnt.lt.0) then
               icnt=-98
               return
            endif
            regs.bregs.dl=char(8)
            regs.bregs.ah=2
            call intdosqq(regs,regs)
          endif
c*****RIGHT ARROW - ECHO THE CHARACTER PREVIOUSLY THERE
          if(ich.eq.77) then
            icnt=icnt+1
            regs.bregs.dl=str(icnt:icnt)
            regs.bregs.ah=2
            call intdosqq(regs,regs)
          endif
c*****PAGE DOWN
          if(ich.eq.81) then
            icnt = -81
            return
          endif
c*****PAGE UP
          if(ich.eq.73) then
            icnt = -73
            return
          endif
        else
          ich=regs.bregs.al
c*****8=BACKSPACE;13=CARRIAGE RETURN;27=ESC
          if(ich.eq.13) then
            if(icnt.eq.0) icnt = -13
            return
          endif
          if(ich.eq.27) go to 20
          if(ich.eq.8) then
            icnt=icnt-1
            if(icnt.lt.0) then
               icnt=-98
               return
            endif
            regs.bregs.dl=char(8)
          else
            icnt=icnt+1
            str(icnt:icnt)=char(ich)
            regs.bregs.dl=str(icnt:icnt)
          endif
          regs.bregs.ah=2
          call intdosqq(regs,regs)
        endif
10    continue
20    if(icnt.eq.0) icnt=-99
      return
      end
c***************************************
      subroutine page(set)
      implicit integer*2 (a-z)
      include 'fgraph.fd'
      character*3 set
      if(set.eq.'set') then
      oldapage=getactivepage()
      oldvpage=getvisualpage()
      newapage=setactivepage(oldapage+1)
      newvpage=setvisualpage(oldvpage+1)
      call clearscreen($GCLEARSCREEN)
      else
      oldapage=setactivepage(oldapage)
      oldvpage=setvisualpage(oldvpage)
      endif
      return
      end
c***************************************
      subroutine edit(ierr)
      implicit integer*2 (a-z)
      parameter(ninp=25,numt=50)
      include 'fgraph.fd'
      record/rccoord/curpos
      integer*4 dum1,bgd,bgd1,bgd2,is,ie,ierr
      logical quest
      character*80 tinp(numt)
      character*80 inp(ninp)
      character*80 string
      character*6 tindex(numt)
      character*1 nul,resp
      logical array
      common/lset/setrow,setcol,linptr,line,inum,inp
      common/lout/tindex,ix,tinp,array
      data bgd/4/fgd/15/,bgd1/1/,fgd1/14/,bgd2/2/,nul/#00/
      ierr = 0
      quest=.false.

* check for mouse, resetting it in the process
5     m1 = 0
      call mouse(m1, m2, m3, m4)

* Quit if mouse wasn't found
      if (m1.eq.0)  then
        print *,  ' Mouse not found'
        ierr=1
        return
      endif

* Initialize menu pointer to first option
      linptr = 1

c*****READ INPUT 23 LINES AT A TIME
      is = 1
      ie=23
c*****DELAY AFTER MOUSE BUTTON PRESSED
c*****SEEMS TO PREVENT ACCIDENTAL PAGE FLIPPING
66    call sleep(.5)
      call input(is,ie,ierr)
      if(ierr.ne.0) go to 300
      nlines=ie
      if(is.eq.1) go to 77
      if(is.eq.2) go to 777
      go to 88
77    write(string,'(a)')'PAGE DOWN for next page or ESC to exit'
      go to 9
777   write(string,'(a)')'PAGE UP for previous page or ESC to exit'
      go to 9
88    write(string,'(2a)') 'press MOUSE or ENTER to return',
     & ' or ESC to exit'
9     dum1= setbkcolor(bgd1)
      call clearscreen($GCLEARSCREEN)
      dum1= setbkcolor(bgd)
      dum2=settextcolor(fgd)
      call settextposition(25,1,curpos)
      call outtext(string(1:42))
      string=' '
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)

c*****RESET MOUSE
      m1 = 0
      call mouse(m1, m2, m3, m4)

* Initialize count of accumulated vertical mouse motion
      motion = 0

* Set flag to update menu first time
      wflag = 1

* Main loop starts here
   10 continue

* Update the menu only when necessary
      if (wflag.eq.1) then
          wflag = 0

        line=0
        do 20 inum = 1,ie
        line=line+1
        call outp
20      continue
        call settextposition(setrow,setcol,curpos)

* End of updating the menu
      end if

c*****CHECK ON KEYBOARD INPUT
      cnt=0
      inum=linptr
      call getstr(inp(inum)(setcol:),cnt)


      if(cnt.ne.0) then
c*****CHECK FOR ENTER KEY AS FIRST CHARACTER
       if(cnt.eq.-13) go to 40
c*****CNT = -99 MEANS ESC AS FIRST TYPED CHARACTER
       if(cnt.eq.-99) go to 999
c*****CNT = -81 MEANS PAGE DOWN(NEXT PAGE)
       if(cnt.eq.-81.and.is.eq.1) then
          is = is + 1
          linptr = 1
          go to 66
       endif
c*****CNT = -73 MEANS PAGE UP(PREVIOUS PAGE)
       if(cnt.eq.-73.and.is.eq.2) then
         is = is -1
         linptr = 1
         go to 66
       endif
c*****CNT = -98 MEANS CURSOR WAS MOVED TO = SIGN OR LEFT OF IT
       if(cnt.eq.-98) go to 30
c*****CNT = -72 MEANS UP ARROW
       if(cnt.eq.-72) then
         if (linptr.gt.1) then
           linptr = linptr - 1
           wflag = 1
         endif
         go to 10
       endif
c*****CNT = -80 MEANS DOWN ARROW
       if(cnt.eq.-80) then
         if (linptr.lt.nlines) then
           linptr = linptr + 1
           wflag = 1
         endif
         go to 10
       endif
c*****CHECK FOR COMMENT OR QUESTION LINE 
       if(inp(inum)(1:1).eq.'#'.or.inp(inum)(1:1).eq.'?'
     & .or.inp(inum)(1:1).eq.nul) go to 30
       line=linptr
       call outp
       call settextposition(setrow,setcol,curpos)
       cnt=0
         ix=ix+1
         tinp(ix)=inp(inum)(setcol:)
         tindex(ix)=inp(inum)(1:6)
       call output('n',ierr)
       if(ierr.ne.0) return
      endif


* Accumulate vertical mouse motion counts
30    m1 = 11
      call mouse(m1, m2, m3, m4)
      motion = motion + m4

* Move up the menu if enough mouse motion
      if (motion.lt.-17) then
          motion = 0
          if (linptr.gt.1) then
              linptr = linptr - 1
              wflag = 1
          end if
      end if

* Move down the menu if enough mouse motion
      if (motion.gt.17) then
          motion = 0
          IF (linptr.lt.nlines) then
              linptr = linptr + 1
              wflag = 1
          end if
      end if

      m1=3
      m2=0
      call mouse(m1, m2, m3, m4)
      if(m2.eq.0) go to 10

40     if(inp(inum)(1:1).eq.'?') then
         isave = is 
         lsave = linptr
         read(inp(inum)(2:3),'(i2)')  is
         is = is + 2
         linptr = 1
         quest=.true.
         go to 66
       endif

      if(quest) then
        quest=.false.
        is = isave
        linptr = lsave
        go to 66
      endif

      go to 10

999   dum1= setbkcolor(bgd2)
      dum2=settextcolor(fgd)
      string=' '

      call settextposition(25,1,curpos)
      write(string,'(a)')'WANT TO SAVE CHANGES TO A DISK FILE ? [y]'
      call outtext(string(1:77))
      call settextposition(25,42,curpos)

      read(*,'(a)') resp
      if(resp.ne.'n'.and.resp.ne.'N') then
c*****SAVE 
          call output('y',ierr)
      endif

300   dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      return
      end
c***************************************
      subroutine sleep(sec)
      integer*2 ihr,imin,isec,i100th
      call gettim(ihr,imin,isec,i100th)
      t1 = ihr*3600. + imin*60. +  isec + i100th/100. 
10    call gettim(ihr,imin,isec,i100th)
      t2 = ihr*3600. + imin*60. +  isec + i100th/100.
      if((t2-t1).lt.sec) go to 10
      return
      end
c***************************************
      subroutine outp
      implicit integer*2 (a-z)
      parameter(ninp=25)
      include 'fgraph.fd'
      integer*4 dum1,bgd,bgd1
      character*80 inp(ninp)
      common/lset/setrow,setcol,linptr,line,inum,inp
      data bgd/7/fgd/0/,bgd1/1/,fgd1/14/

* Write menu, highlighted if selected
      if (linptr.eq.line) then
      dum1= setbkcolor(bgd)
      dum2=settextcolor(fgd)
      else
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      end if
      call placel
      return
      end
c***************************************
      subroutine placel
      implicit integer*2(a-z)
      include 'fgraph.fd'
      record/rccoord/curpos
      parameter(ninp=25)
      character*80 inp(ninp)
      common/lset/setrow,setcol,linptr,line,inum,inp
      character*3 str
c*****STR CONTAINS HEX CODES FOR =  ? and NUL
      data str(1:1)/#3d/,str(2:2)/#3f/,str(3:3)/#00/
      call settextposition(line,1,curpos)
      call outtext(inp(inum))
      if(linptr.eq.line) then
      setrow= line
      setcol = scan(inp(inum),str) + 1
      endif
      return
      end
c***************************************
      subroutine output(resp,ierr)
c
c     namelist output
c
      parameter(numt=50)
      include 'fgraph.fd'
      record/rccoord/curpos
      character*1 resp
      integer*4 dum1,bgd1
      integer*2 ix
      character*6 tindex(numt)
      character*80 tinp(numt)
      logical array
      common/lout/tindex,ix,tinp,array
      character str*80,eq*4
      character*40 sfile
      logical save
      data bgd1/1/,fgd1/14/,jcmd/9/,eq/' = '/
c
 
      ierr=0

      if(resp.eq.'y') then
 
      dum1= setbkcolor(bgd1)
      dum2=settextcolor(fgd1)
      call clearscreen($GCLEARSCREEN)

c*****SAVE
        if(ix.eq.0) then
          call settextposition(10,1,curpos)
          call outtext('NO VALUES WERE INPUT')
          call settextposition(11,1,curpos)
          ierr=1
          return
        endif

        call settextposition(10,1,curpos)
        write(*,'(a\)') ' ENTER FILE NAME FOR SAVING VALUES: '
        read(*,'(a)') sfile
        open(unit=jcmd,file=sfile,form='formatted',
     1  status='unknown',err=992,iostat=icheck)
        save = .true.
      else
c*****DON'T SAVE
        save = .false.
      endif
      go to 50

c*****FILE ERRORS
992   write(str,'(a,i5,2a)') 'error no. ',icheck,
     & ' on opening file: ',sfile
      go to 995
993   write(str,'(a,i5,2a)') 'error no. ',icheck,
     & ' on writing file: ',sfile
      close (jcmd)
995   call settextposition(11,1,curpos)
      call outtext(str(1:75))
      sfile=' '
      call settextposition(12,1,curpos)
      ierr=1
      return

c*****PUT VALUES IN A DISK FILE OR ARRAY TINP
c*****TO BE READ INTO MEMORY BY SUBROUTINE NAMEL
c*****FIRST CHECK FOR DUPLICATE LINE CHANGES AND USE ONLY LAST ONE
50    it = 0
      do 52 i=ix,1,-1
      do 51 j=i-1,1,-1
      if(tindex(i).eq.tindex(j)) tindex(j)=' '
51    continue
52    continue
      if(save) then
        write(jcmd,'(a)',err=993,iostat=icheck) ' &value'
      endif
     
      do 100 i=1,ix
      if(tindex(i).eq.' ') go to 100
c******VALUES ON ONE LINE
       if(save) then
         write(jcmd,'(1x,2a,a58)',err=993,iostat=icheck)
     &        tindex(i),eq,tinp(i)
       else
         it = it + 1
         write(str,'(a)') tinp(i)(1:70)
         write(tinp(it),'(3a)') tindex(i),eq,str(1:70)
c         print *,'output ',it,i,tindex(i),' gap ',tinp(it)
c         pause
       endif
100   continue
      if(save) then
        write(jcmd,'(a)',err=993,iostat=icheck) ' &'
      else
         it = it + 1
         tinp(it) = ' &'
      endif
      if(save) then
        rewind jcmd
        call namemc(jcmd,ierr)
        close(jcmd)
      else
        call namemc(-99,ierr)
      endif
      return
      end
c************************************************************************
      subroutine inptxt(is,ie)
      character*80 inp(25)
      integer*2 setrow,setcol,line,linptr,inum
      common/lset/setrow,setcol,linptr,line,inum,inp
      character*1 ap,nul
      character*56 ifile,ofile,title,ifmt
      common/fmt/ifmt,izchan,idirx,maxrec
      common/gparm/mxc,mxr,nc,nr,nsec,ntier,lapovr,
     & nim,epsm,dv,slope,gmin,gmax
      common/nlparm/ifile,ofile,title,xo,yo,del,nco,
     & nro,npmin,radius,hull,whole,region,iproj,cmdeg,bldeg
      data ap/''''/,iptr/0/,nul/#00/

      go to (1,2),is
1     ii=iptr+1
      write(inp(ii),'(3a)') nul,'THE FOLLOWING SEVEN PARAMETERS',
     & ' ARE REQUIRED' 
      ii=ii+1
      write(inp(ii),'(a)')nul 
      ii=ii+1
      write(inp(ii),'(4a)') 'ifile (INPUT FILE)=',ap,
     & ifile,ap
      ii=ii+1
      write(inp(ii),'(4a)') 'ofile (OUTPUT FILE)=',ap,
     & ofile,ap
      ii=ii+1
      write(inp(ii),'(2a,t60,a,g15.6)') 'xo    (X(E-W) COORDINATE',
     & ' OF THE LOWER-LEFT CORNER OF GRID)','=',xo
      ii=ii+1
      write(inp(ii),'(2a,t60,a,g15.6)') 'yo    (Y(N-S) COORDINATE',
     & ' OF THE LOWER-LEFT CORNER OF GRID)','=',yo
      ii=ii+1
      write(inp(ii),'(2a,t60,a,g15.6)') 'del   (GRID INCREMENT',
     & ' THAT APPLIES TO BOTH X & Y DIRECTIONS)','=',del
      ii=ii+1
      write(inp(ii),'(2a,t60,a,i5)') 'nc    (NUMBER OF OUTPUT',
     & ' COLUMNS; X DIRECTION)','=',nc
      ii=ii+1
      write(inp(ii),'(2a,t60,a,i5)') 'nr    (NUMBER OF OUTPUT ROWS',
     & '; Y DIRECTION)','=',nr
      ii=ii+1
      write(inp(ii),'(a)')nul 
      ii=ii+1
      write(inp(ii),'(3a)') nul,'     ADDITIONAL OPTIONAL PARAMETERS',
     & ' ARE AVAILABLE ON THE NEXT PAGE' 
      ie=ii
      return
2     ii = iptr + 1
      write(inp(ii),'(4a)') 'title (56 CHARACTERS)=',ap,
     & title,ap
      ii=ii+1
      write(inp(ii),'(2a,t60,a,f6.0)') 'radius (DISTANCE TO',
     & ' DETERMINE WHERE NO-DATA FLAGS ARE PUT)', '=',radius
      ii=ii+1
      write(inp(ii),'(2a,t60,a,f6.0)') 'slope (DISTANCE WEIGHTING',
     & ' FACTOR)', '=',slope
      ii=ii+1
      write(inp(ii),'(a,t60,a,f6.0)')  'hull  (FILLS NO DATA AREAS)',
     & '=',hull
      ii=ii+1
      write(inp(ii),'(2a,t60,a,g15.6)')'gmin  (MINIMUM DATA VALUE',
     & ' TO USE FOR GRIDDING)','=',gmin
      ii=ii+1
      write(inp(ii),'(2a,t60,a,g15.6)') 'gmax  (MAXIMUM DATA VALUE',
     & ' TO USE FOR GRIDDING)','=',gmax
      ii=ii+1
      write(inp(ii),'(2a,t60,a,i5)') 'npmin (NO. OF VALUES',
     & ' < radius BEFORE A NO-DATA FLAG IS PUT)','=',npmin
      ii=ii+1
      write(inp(ii),'(2a,t60,a,f6.0)') 'epsm  (MAXIMUM CHANGE PER',
     & ' INTERATION)','=',epsm
      ii=ii+1
      write(inp(ii),'(2a,t60,a,i5)') 'nim   (MAXIMUM NUMBER OF',
     & ' INTERATIONS PER BLOCK)','=',nim
      ii=ii+1
      write(inp(ii),'(2a,t60,a,i5)') 'idirx (SET TO  NON-ZERO',
     & ' WHEN USING POSITIVE WEST LONG.)','=',idirx
      ii=ii+1
      write(inp(ii),'(2a,t60,a,f6.0)') 'region (SET NON-ZERO TO SAVE',
     & ' REGIONAL GRID AS REGIONAL.TMP)','=',region
      ii=ii+1
      write(inp(ii),'(2a,t60,a,f6.0)') 'whole (SET NON-ZERO TO SAVE',
     & ' GRID BEFORE radius USED)','=',whole
      ie=ii
      return
      end
c************************************************************************
c                            Appendix B
c                    USGS Open-File Report 81-1224
c
c                minimum curvature gridding routine
c************************************************************************
c  This program generates a 2-dimensional grid, equally
c  incremented in x and y, from randomly placed data points.
c  The algorithm (Briggs) produces a smooth grid by iteratively
c  solving a set of difference equations which minimize the total
c  2nd horizontal derivative and attempt to honor input data.
c  (ref: I.C. Briggs, 1974, Geophysics, v 39, no 1)
c
c  namelist parameters:
c
c       nonoptional parameters
c  xo      x coordinate of lower left corner of grid
c  yo      y coordinate of lower left corner
c  del     x and y increment (must be positive).  xo, yo, del, and the
c          input x and y data points are in the same coordinate system.
c  nc      number of columns
c  nr      number of rows ( nc*nr < 1.3e6 )
c
c       parameters which can be interactive   
c  ifile   input filename containing xyz data records.  see subroutine
c          'openf' for complete description of recognized structures.
c  ofile   output grid filename, consisting of a header record 
c          and gridded z values in row records.
c  title   56 character title of output grid entered in the 
c          header record. (default is blank)
c  ifmt    input format: present if input is ascii. default is blanks
c          indicating either binary input data or free field ascii.
c  izchan  selection of z channel from an input record of the form:
c          station_id(8char), x, y, channel(1), channel(2)...channel(6).
c          Internally, izchan is set to izchan+2, and channels 1 and 2
c          are assumed to be x and y, as this program was originally written.
c          Default = zero, xyz data.
c
c       optional parameters
c  radius  in horizontal data units, grid points with no data inside
c          this radius have a 'no data' value inserted (dval).
c          default = 0, ie. infinity - a completely defined grid.
c  hull    apply hull algorithm to eliminate dval islands ( no-data
c          areas completely surrounded by data ), note radius must be
c          set.  hull is a switch: 0=off, 1=on, default=off.
c  slope   a distance weighting function to decrease aliasing by
c          combining all data in a grid cell to produce one data 
c          value. (default = 0, only the nearest data effect grid nodes).
c          slope contols a bell function: weight = 1 / [ dist^2 + 1/slope ]
c          where distance is in units of grid cells. 
c  gmin/gmax  minimum and maximum grid values allowed, in z data units.
c          resistivity data might have gmin = 0 for example.
c          defaults: gmin=-1.e37,  gmax=1.e37
c
c  npmin   number of data points within 'radius' distance
c          before grid point considered valid. (default = 1)
c  epsm    in z data units, iteration cutoff (default = 0).
c  nim     maximum iterations per block (default = 20).
c  idirx   set to 1 when x coordinates decrease with increase in column,
c          postive west longitude for instance. (default = 0)
c  region  set .ne. 0. to save regional grid
c  whole   set .ne. 0. to save unradiused grid for subsequent operators
c          which require completely defined grid.
c
c    The next three parameters allow subsequent programs to recover 
c    map projection information.  Minc does not use these numbers.
c
c  iproj   site dependent map projection number
c  cm      central meridian in degrees
c  baslat  base latitude in degrees
c
c
c               algorithm description
c
c  Minc breaks the grid area into blocks containing no more than
c  5000 points.  For each block: a temporary binary file containing the
c  input data is read, an initial grid is interpolated using
c  one-dimensional interpolation to fill holes, data points are assigned
c  to grid points, and iteration using minimum curvature difference
c  equations attempts to honor the data points.
c  Continuity between blocks is maintained by initializing the next block
c  with interpolated values from a coarser regional surface and with a 
c  given amount of overlap from preceeding blocks.
c
c                    Mike Webring
c                 U.S. Geological Survey
c                 Geophysics branch, mail stop 964
c                 PO box 25046
c                 Denver Federal Center
c                 Denver, Colorado 80225
c
c       main structure Oct. 1978
c       VAX version March 83
c
c       subroutine ASSIGN,CELL3,ZGINT 2/84
c smooth interpolation of the regional as a starting solution for 
c the final grid.
c
c       EXPRESS file input 2/87
c along with hull algorithm and automatic deciphering of the input data
c file structure.  many routines converted to fortran 77. 
c
c
c   Disclaimer:
c   Although this program has been used by the U.S. Geological Survey,
c   no warranty, expressed or implied, is made by the Survey as to the
c   accuracy and functioning of the program and related program
c   material nor shall the fact of distribution constitute any such 
c   warranty, and no responsibility is assumed by the Survey in 
c   connection therewith.
c
c************************************************************************

       subroutine mincs(ierr)
	parameter       ( maxchn=30, nwork=40000, npidg=375, npb=100 )
	common /iodev/  igrid, iregn, imask, idata, irand, isubg
	common /nlparm/ name1(14), name2(14), name3(14), xo, yo, del,
     1               nco, nro, npmin, radius, hull, whole, region,
     1               iproj, cmerid, baslat
	common /fmt/    ifmta(14), izchan, idirx, maxrec
	common /datai/  xyz(3*npb)
	common /contn1/ zl(npidg), zb(npidg)
	common /rgnl/   ifill(4), mdel, rg(5000)
	common /gparm/  mxc, mxr, nc, nr, nsec, ntier, lapovr, 
     1               nim, epsm, dval, slope, gmin, gmax
	common /qparm/  ihwind, mxcq, mxrq, maxq, maxg
	common /assem/  ntot(2)
	common /array/  work(nwork)

	dimension iwork(nwork)
	dimension       za(maxchn), idim(3), xloc(4)
	character*56    ifile, ofile, title, ifmt, blank
	character*8      p
	logical         binary, packed,exists,opn
	equivalence(iwork(1),work(1))
	equivalence     ( ifile, name1 ), ( ofile, name2 ),
     1               ( title, name3 ), ( ifmt,  ifmta )

       ierr = 0
	blank  = ' '
	p   = 'min-curv'
	work(nwork) = 0.
	zb(npidg)   = 0.
	npb3        = 3*npb
	xyz(npb3)   = 0.
	mdel    = 0
	maxg    = 0
	ntot(1) = 0
	ntot(2) = 0

c  extend grid so general equation can be used 
c  at the edge of the desired grid.
	mxc = nco + 4
	mxr = nro + 4

c  open input data file and derive structure.

       if(ifile.ne.'tr.dat') then
         call openf( idata, ifile, izchan, iwork, maxchn, 
     1           binary, packed, nword, ndmax, ntotal, ierr )
         if(ierr.ne.0) go to 9999
c      if ( .not. binary .and. ifmt(1:4) .eq. '    ' ) then
c        print *, ' ascii input data, free field formatting'
c      endif
       else
       inquire(file='tr.dat',exist=exists)
       if(exists) then
         inquire(file='tr.dat',opened=opn)
         if(.not.opn) then
           open(idata,file='tr.dat',status='old',err=900,iostat=icheck)
           read(idata,'(14x,a)',err=900,iostat=icheck) 
           read(idata,*,err=900,iostat=icheck) dum,dum,dum,dum
           read(idata,*,err=900,iostat=icheck) dum,idum,idum
           go to 3
900        print *,'error no. ',icheck,' with file tr.dat'
           ierr = 1
           go to 9999
         endif
       else
         print *,'file tr.dat does not exist'
         print *,'run program again with augumentation'
         ierr = 1
         go to 9999
       endif
c******FOLLOWING VARIABLES LIFTED FROM openf TO REFLECT XYZ ASCII INPUT
         nword  = 0
	  ndmax  = 0
         ntotal = 0
         binary = .false.
         packed = .false.
         izchan = 0
       endif


c  check work space required for masking

3      ihwind = int( abs( radius / del ) + 0.5 )
	nreq   = (nco + 2*ihwind) * ( 2*ihwind + 3 )
	if ( nreq .gt. nwork ) then
	  print *,' radius too large for available work array.' 
	  b        = 2.0 * nco + 6.0
	  c        = 3.0 * nco - nwork
	  ihwind   = int( (sqrt( b*b - 16.0*c ) - b) / 8.0 )
	  iwind    = 2 * ihwind + 1
	  radmax   = del * float(ihwind)
	  print 5, nco, radmax
5	  format('  maximum radius for a grid with ',i6,' columns is ',
     1        /, 1pe14.5, ' data units' )
	  print *, ' continuing with radius of ', ihwind, 
     1         ' grid cells'
	endif
c  open output grid.
c  VAX allows fixed length binary files to be read either sequential
c  or direct.  OFILE is created sequentially.
	length = nco + 1
	if ( length .lt. 30 ) length = 30
	if ( ofile .eq. blank ) then
	  ofile = 'trminc.grd'
	  print *, ' output grid filename is "trminc.grd" '
	endif
c  open output file depending on value of ihwind
      if(ihwind.gt.0) then
      open(igrid,status='scratch',form='unformatted')
      else
      open(igrid,file=ofile,status='unknown',
     & form='unformatted')
       endif

c  adjust coordinates and write grid header record.
       nc = 0
       nr = 0
	delx = del
	xout = xo
	if ( idirx .ne. 0 ) then
	  delx  = -del
	  idirx = -1
	  xo    = -xo
	endif
	nz = 1
       write( igrid ) title,p, nco,nro,nz, xout,delx,yo,del, 
     1         iproj, cmerid, baslat

c  setup masking grid

	if ( ihwind .gt. 0 ) then
	  iwind  =   1 + 2   * ihwind
	  mxcq   = nco + 2   * ihwind
	  mxrq   = nro + 2   * ihwind
	  xo2    = xo  - del * ihwind 
	  yo2    = yo  - del * ihwind 
	  length = mxcq + 1
	  if ( length .lt. 30 ) length = 30
	open(unit=imask,status='scratch',
     1       form='unformatted', recl=length*4, access='direct' )
	  write( imask, rec=1 ) title,p, mxcq,mxrq,nz, xo2,del,
     1        yo2, del, iproj, cmerid, baslat
	endif

c  setup grid partition parameters

	call prtish( nwork, npidg, mxc, mxr, 
     1            nc, nr, nsec, ntier, lapovr )
	print  10, nc, nr, nsec, ntier
10	format(/,' partition size:          ncol=',i3,
     & ',         nrow =',i3,/,
     1 ' arrangement:  column sections=',i3,',  row sections=',i3,/)
	if ( nsec .eq. 0 ) then
          print *,' number of sections is zero'
          ierr=1
          go to 9999
       endif
	ng   = nc*nr
	nblk = nsec*ntier
	if ( nblk .gt. npidg ) then
          print *, 'argh...grid too large'
          ierr=1
          go to 9999
       endif
	n    = int( float(nwork) / float(3*nblk) )
	if ( n .gt. npb )  n = npb
	npb3 = 3 * n
	if ( nblk .eq. 1 ) lapovr = 1

c  read data, prepare random access file.
c  record length in tmp is 2+npb3 words, the number of records depends
c  on input data file.  

	length = 2 + npb3
	open ( unit=irand, access='direct',
     1     form='unformatted', status='scratch',
     1     recl=length*4 )
	xo2 = xo - 2.*del
	yo2 = yo - 2.*del

	if ( packed ) then
c  read express file
	  call expdat( idata, irand, del, xo2, yo2, npb3, 
     1              rg, nword, ndmax, ntotal, ierr )
         if(ierr.ne.0) go to 9999
	  else
c  single data point per record.
	  lenz = nword - 2
	  if ( lenz .le. 0 ) lenz = 1
	  call seqdat( idata, irand, del, xo2, yo2, npb3, 
     1              binary, za, lenz, ierr )
         if(ierr.ne.0) go to 9999
	endif

       inquire(file='tr.dat',exist=exists)
       if(exists) then
         inquire(file='tr.dat',opened=opn)
         if(opn) then
           close(idata,status='delete')
         else
           open(16,file='tr.dat',status='old')
           close(16,status='delete')
         endif
       endif
	close( idata )

c  record length in tmp2 is nc*nr words
c  maximum number of records is 2*nsec
	open (unit=isubg, access='direct', 
     1     form='unformatted', status='scratch',
     1      recl=nc*nr*4)

c  regional grid
	if ( region .ne. 0.0 ) open( unit=iregn, file='regional.tmp',
     1                      status='unknown', form='unformatted' )
	call rejonl( region, xo2, yo2, del, xyz, npb3, ierr )
         if(ierr.ne.0) go to 9999

c  produce grid at specified interval.
	call pcontl( nco, nro, ng, work, npb3, xyz, ierr )
         if(ierr.ne.0) go to 9999

	close( isubg )
	close( irand )
	close( iregn )

c  trim completed grid to data coverage

	if ( ihwind .eq. 0  ) go to 999

       if ( whole .ne. 0.0 ) then
c  copy completely defined grid into file whole.tmp
	  ifull = 17
	  open( unit=ifull, file='whole.tmp', status='unknown', 
     1       form='unformatted' )
	  rewind( igrid )
	  read ( igrid ) title, p, idim,xloc, iproj,cmerid,baslat
	  write( ifull ) title, p, idim,xloc, iproj,cmerid,baslat
	  do 20 j = 1, idim(2)
	    call rowio( idim(1), iwork, 1, igrid, ifull, ie )
   20     continue
	  close( ifull )
       endif

c  apply radius operator to masking grid
	n1 = 1  + mxcq * iwind
	n2 = n1 + mxcq
	call radmsk( npmin, imask, mxcq, iwind, 
     1            iwork, iwork(n1), iwork(n2), ierr )
         if(ierr.ne.0) go to 9999

	if ( hull .ne. 0 ) then
c  fill no-data polygons that intersect the grid border with flag=2.
	  call outside( imask, ihwind, iwork, nwork, ierr )
         if(ierr.ne.0) go to 9999
	  iflag = 2
	  else
c  all no-data polygons
	  iflag = 0
	endif

c  use masking template to insert dval's
c      close( igrid )
c      open( igrid, file=ofile, status='old', form='unformatted',
c     1     access='direct' )
        rewind igrid
	  read ( igrid ) title, p, idim,xloc, iproj,cmerid,baslat
	open(unit=igrid+10,file=ofile,status='unknown',form='unformatted')
	  write ( igrid+10 ) title, p, idim,xloc, iproj,cmerid,baslat
	n1     = mxcq + 1
	istart = ihwind + 1
	call templat( istart, istart, iflag, dval, 
     1             imask,     iwork, mxcq, mxrq,
     1             igrid, work(n1),  nco,  nro, ierr )
       if(ierr.ne.0) go to 9999
	close( igrid )
	close(igrid+10)
	close( imask )

999	continue
       nc=nco
       nr=nro
       xo=xout
       if(izchan.ne.0) izchan=izchan-2
       if(slope.ne.0.) slope = 1.0/slope
9999   if(ierr.ne.0) return
       print *,' '
       pause 'press ENTER to return to menu'
	return
	end
c************************************************************************
	subroutine input(is,ie,ierr)

c  initialize minc controlling parameters
c  The namelist is non-standard Fortran-77 and provides freeform input
c  suitable for program command files where only the necessary parameters
c  are included.
c  An example command file would be:
c     &parms
c     ifile='test.dat', ofile=' test.grd'
c     xo=10., yo=0., del=.25, nc=100, nr=75
c     radius=1.25, hull=1
c     &

	parameter       ( maxchn=30, nwork=40000, npidg=375, npb=100 )
	common /iodev/  igrid, iregn, imask, idata, irand, isubg
	common /nlparm/ name1(14), name2(14), name3(14), xo, yo, del,
     1               nco, nro, npmin, radius, hull, whole, region,
     1               iproj, cmdeg, bldeg
	common /fmt/    ifmta(14), izchan, idirx, maxrec
	common /gparm/  mxc, mxr, nc, nr, nsec, ntier, lapovr, 
     1               nim, epsm, dval, slope, gmin, gmax
	common /gparm2/ cm(3),baslat(3)
	common /qparm/  ihwind, mxcq, mxrq, maxq, maxg
	common /assem/  ntot(2)

	
c      namelist /parms/ title, nc, nr, xo, yo, del, idirx,
c    1c    c    c     ifile, ofile, ifmt, izchan,
c    1c    c    c     radius, npmin, hull, slope,
c    1c    c    c     gmin, gmax, iproj, cm, baslat,
c    1c    c    c     nim, epsm, whole, region

	character*56    ifile, ofile, title, ifmt, blank, cfile
	equivalence     (ifile, name1),    ( ofile, name2), 
     1               ( title, name3 ),  ( ifmt, ifmta )
c  dval (default grid value) is the no-data flag, approximate value = 1.0e38.
       logical first,exists
      parameter(numt=50)
      character*80 tinp(numt)
      character*6 tindex(numt)
      character*1 ap
      integer*2 ix
      logical array
      common/lout/tindex,ix,tinp,array
	data  dv/1.0e38/,first/.true./,ap/''''/

       ierr = 0
	icmd   = 9
	blank  = ' '

c  iodev common
	igrid = 10
	iregn = 11
	imask = 12
	idata = 13
	irand = 14
	isubg = 15

       if(first.or.is.eq.50) then
       ix = 0
c  nlparm common
       ifile  = ' '
       ofile  = ' '
	title  = ' '
	del    = 0.
	npmin  = 1
	radius = 0.0
	hull   = 0.0
	region = 0.0
	whole  = 0.0
	iproj  = 0
c      do 10 i = 1, 3
c      cm(i)      = 999.0
c      baslat(i)  = 999.0
c   10   continue

c  fmt common
	ifmt   = ' '
	izchan = 0
	idirx  = 0

c  gparm common
	nc     = 0
	nr     = 0
	lapovr = 10
	nim    = 20
	epsm   = 0.0
	dval   = dv
	slope  = 0.0
	gmin   = -1.0e37
	gmax   =  1.0e37
      inquire(file='tr.dat',exist=exists)
      if(exists.and.first) then
        open(idata,file='tr.dat',status='old',err=900,iostat=icheck)
        read(idata,'(14x,a)',err=900,iostat=icheck) ofile
        read(idata,*,err=900,iostat=icheck) xo,dum,yo,dum
        read(idata,*,err=900,iostat=icheck) del,nc,nr
        ifile  = 'tr.dat '
         ix=ix+1
         tinp(ix)=ap//ifile//ap
         tindex(ix)='ifile '
         ix=ix+1
         tinp(ix)=ap//ofile//ap
         tindex(ix)='ofile '
         ix=ix+1
         write(tinp(ix),'(g15.6)') xo
         tindex(ix)='xo    '
         ix=ix+1
         write(tinp(ix),'(g15.6)') yo
         tindex(ix)='yo    '
         ix=ix+1
         write(tinp(ix),'(g15.6)') del
         tindex(ix)='del   '
         ix=ix+1
         write(tinp(ix),'(i5)') nc
         tindex(ix)='nc    '
         ix=ix+1
         write(tinp(ix),'(i5)') nr
         tindex(ix)='nr    '
      endif
      endif
       first=.false.
       if(is.eq.50) go to 1
       if(is.eq.10) go to 3
       call inptxt(is,ie)
       return

c  read command file.

1     write(*,'(a\)') ' ENTER COMMAND FILE NAME: '
      read(*,'(a)') cfile
      open(unit=icmd,file=cfile,form='formatted',
     1 status='old',err=2,iostat=icheck)
      call namemc(icmd,ierr)
      close (icmd)
      return
2     call fcheck(cfile,icheck,ierr)
      if(ierr.eq.0) go to 1
      if(ierr.ne.0) then
        close (icmd)
        return
      endif

ccc
ccc izchan modified to accommodate internal convention
3     if(izchan.gt.0)izchan=izchan+2
ccc

c  check and initialize basic parameters

	if ( nc .le. 4  .or.  nr .le. 4 ) then
          print *, ' parameter NC or NR is less than 5'
          ierr=1  
          return
       endif
	if ( del    .le. 0.0 )   then
          print *,' parameter DEL is less than or equal zero'
          ierr=1  
          return
       endif
	if ( lapovr .lt. 4   )   lapovr =  4
	if ( slope  .ne. 0.0 )   slope  =  1.0 / abs( slope )
	if ( radius .lt. 0.0 )   radius =  0.0
	if ( hull   .ne. 0.0 ) then
	  if ( radius .le. 0.0 ) 
     1      print *, ' %%radius must be set for hull to operate'
	  if ( radius .gt. 0.0  .and.  radius .le. del ) 
     1      radius = 1.1 * del
	endif
	nco = nc
	nro = nr

	if ( ifile .eq. blank ) then
         write(*,'(a\)') ' enter xyz data filename :'
	  read 37, ifile
37       format(a)
         ix=ix+1
         tinp(ix)=ap//ifile//ap
         tindex(ix)='ifile '
	endif
c       write(*,'(a\)') ' enter xyz data format or <car. ret.>:'
c       read 37, ifmt
c       if ( ifmt(1:4) .eq. '    ' ) then
c       write(*,'(a)') ' channel number to be gridded '
c       write(*,'(a\)') ' (0 for xyz files, 1-6 for post files) :'
c         read *, izchan
ccc
ccc  izchan modified to accommodate internal convention
c        if(izchan.ne.0)izchan=izchan+2
ccc
c        if ( izchan .gt. maxchn  .or. izchan .lt. 0 ) then
c           write(*,'(a,i3)') ' 0 <= izchan <= ', maxchn
c           ierr=1  
c           return
c         endif
c        endif

	if ( ofile .eq. blank ) then
         write(*,'(a\)') ' enter output grid filename :'
	  read 37, ofile
         ix=ix+1
         tinp(ix)=ap//ofile//ap
         tindex(ix)='ofile '
       endif

c       write(*,'(a\)') ' enter title :'
c       read 37, title

c      call dms2d( cm, cmdeg )
c      call dms2d( baslat, bldeg )

	return
900    print *,'error no. ',icheck,' with file tr.dat'
       ierr = 1
       return
	end
c***************************************
      subroutine fcheck(fname,icheck,ierr)
      character*(*) fname
      character*1 resp
      logical exists
      inquire(file=fname,exist=exists)
        if(.not.exists) then
          print *,'can''t find file ',fname
          fname=' '
          write(*,'(a\)') ' TRY AGAIN ? [y]: '
          read(*,'(a)') resp
          if(resp.eq.'n'.or.resp.eq.'N') then
            ierr=1
            return
          else
            ierr=0
            return
          endif
        endif
      print *,'error no. ',icheck,' with file: ',fname
      fname=' '
      ierr=1
      return
      end
c*****************************************
	subroutine dms2d( dms, d )
	dimension dms(3)
c  take sign from most significant element then convert to decimal degree.
	signd = 1.
	if ( dms(1) .lt. 0.0 ) signd = -1.
	if ( dms(1) .eq. 0.0  .and.  dms(2) .lt. 0.0 ) signd = -1.
	if ( dms(1) .eq. 0.0  .and.  dms(2) .eq. 0.0  .and.
     1    dms(3) .lt. 0.0 ) signd = -1.
	d = sign( dms(1), signd ) + sign( dms(2), signd ) / 60. +
     1   sign( dms(3), signd ) / 3600.
	return
	end
c*****************************************
	subroutine openf( iswt, ifile, izchan, ibuff, maxchn,
     1                 binary, packed, nword, ndmax, ntotal, ierr )
c
	parameter  ( maxbuf=1024 )
	dimension  ibuff(maxbuf)
	character  ifile*(*)
       character*1 ifmtv
	logical    binary, packed
       ierr = 0
       nword  = 0
	ndmax  = 0
	ntotal = 0
	binary = .false.
	packed = .false.

3      call opnsta(iswt,ifile,ifmtv,ierr,*5)
       if(ierr.ne.0) return
       go to 7
5      write(*,'(a\)') ' enter input xyz file name:'
       read(*,'(a)') ifile
       go to 3
c  binary or ascii input ?
7      if(ifmtv.eq.' ') binary = .true. 

999	return
	end
c***********************************************************************
      subroutine opnsta(ista,ifile2,ifmtv,ierr,*)
      character*(*) ifile2,ifmtv
      character*1 resp
      logical exists,opn
      ierr=0
      inquire(file=ifile2,exist=exists)
        if(.not.exists) then
          print *,'can''t find file ',ifile2
          ifile2=' '
          write(*,'(a\)') ' TRY AGAIN ? [y]: '
          read(*,'(a)') resp
          if(resp.eq.'n'.or.resp.eq.'N') then
            ierr=1
            return 
          else
            return 1
          endif
        endif
        inquire(unit=ista,opened=opn)
          if(opn) close(ista)
        open(unit=ista,file=ifile2,
     &  status='old',err=442,iostat=icheck)
        read(ista,*,err=442,end=442,iostat=icheck) x,y,z
        rewind ista
        ifmtv(1:1) = '*'
        return
442     close(ista)
        open(unit=ista,file=ifile2,form='unformatted',
     &  status='old',err=443,iostat=icheck)
        read(ista,err=444,end=444,iostat=icheck) x,y,z
        nword = lenrec(ista,'u')
        if(nword.ne.3.and.nword.ne.-1) go to 444
        ifmtv = ' '
        return
443    print *,'error no. ',icheck,' with file: ',ifile2
       go to 445
444    print *,'file is not an xyz file'
445    ifile2=' '
       close(ista)
       ierr = 1
       return 
       end
c***********************************************************************
	function lenrec( ista,rtype )
c  brute force recovery of record length
       parameter(nwork=50)
	real irec(nwork)
       character*1 rtype

	istep = 10
	rewind ista

	do 5 max = istep, nwork, istep
         if(rtype.eq.'u') then
           read( ista, err=10 ) ( irec(i), i=1, max )
         else
	    read( ista,*, err=10 ) ( irec(i), i=1, max )
         endif
	  rewind ista
5       continue
	lenrec = -1
	go to 99

10	rewind ista
	do 30 j = max, max - istep, -1
         if(rtype.eq.'u') then
	    read( ista, err=20 ) ( irec(i), i=1, j )
         else
	    read( ista,*, err=20 ) ( irec(i), i=1, j )
         endif
	  lenrec=j
	  go to 99
20	  rewind ista
30      continue

99	rewind ista
	return
	end
c******************************************************************************
	subroutine prtish(nwrk,nsave,mc,mr,nc,nr,nsec,ntier,lap)
c  optimize subdivision of the main grid.
c  maximize block size subject to the side ratio less than 3:1
	dimension nc1(3),nr1(3)
c
	nmax=nwrk/8
	is=0
	js=0
	width=aint(sqrt(float(nmax))-float(lap))
	nsec1=int(float(mc)/width+.5)-1
	ntier1=int(float(mr)/width+.5)-1
	if(nsec1.lt.1) nsec1=1
	if(ntier1.lt.1) ntier1=1
	n=nsec1
	do 1 i=1,3
	nc1(i)=int(float(mc-lap)/float(n)+.9999)+lap
1	n=n+1
	n=ntier1
	do 2 i=1,3
	nr1(i)=int(float(mr-lap)/float(n)+.9999)+lap
2	n=n+1
	nb=1000
	do 3 i=1,3
	do 3 j=1,3
	ratio=float(nc1(i))/float(nr1(j))
c  ratio limits match zl&zb array size
	if(ratio.gt.3.0 .or. ratio.lt..3333) go to 3
	if(nc1(i)*nr1(j).gt.nmax) go to 3
c  3 cols and rows are saved to provide continuity between blocks
	if(nc1(i)*3. .gt. nsave) go to 3
	if(nr1(j)*3. .gt. nsave) go to 3
	nblk=(nsec1+i-1)*(ntier1+j-1)
	if(nblk.ge.nb) go to 3
	nb=nblk
	js=j
	is=i
3	continue
	if(is.eq.0) go to 9
	nc=nc1(is)
	nr=nr1(js)
	nsec=nsec1+is-1
	ntier=ntier1+js-1
	return
9	print  8
8	format(' problem with partition')
	nsec=0
	return
	end
c******************************************************************************
	subroutine seqdat( iswt, jswt, del, xo, yo, npb3, 
     1                  binary, za, nz, ierr )

c  Input single record xyz triples.
c  there are nsec*ntier pigeon holes, each of
c  which contains all data necessary for iterating
c  a subgrid.  because of overlap on left and
c  bottom sides one data point can appear in several holes.

	common /array/  wrk(40000)
	common /fmt/    ifmt(14), izchan, idirx, maxrec
	common /gparm/  mxc, mxr, nc, nr, nsec, ntier, lap, nim, epsm,
     1               dval, slope, gmin, gmax
	common /contn1/ loca(375), ioff(375)

	dimension       za(nz)
	character*8     id
	character*4     tstfmt
	equivalence     ( ifmt(1), tstfmt )
	logical         binary

       ierr = 0

	if ( binary ) then
	  if ( izchan .eq. 0 ) then
	    itype = 1
	    else 
	    itype  = 2
	    ixchan = 1
	    iychan = 2
	  endif
	  else
	  if ( tstfmt .ne. '    ' ) then
	    itype = 3
	    else
	    itype = 4
	  endif
	endif

	nxyz = npb3
	jrec = 0
	ic   = 0
	nblk = nsec * ntier
	nchk = npb3 - 3
	npb  = npb3 / 3
	err  = 1.e-2
	fudg = del * err
	xmax = float( mxc - 1 ) * del + xo - fudg
	ymax = float( mxr - 1 ) * del + yo - fudg
	xmin = xo + fudg
	ymin = yo + fudg

c  dimension of pigeon hole minus overlap
	cl   = float( nc - lap )
	rl   = float( nr - lap )
	dux  = 1. / cl
	duy  = 1. / rl
	fx   = ( float( lap ) - 0.1 ) * dux
	fy   = ( float( lap ) - 0.1 ) * duy

c  ensures smallest index (igx,y) in 'bwts' is 1
	err2 = 1.0 + 0.5 * err
	xlmt = 1.0 - err2 * dux
	ylmt = 1.0 - err2 * duy
	rdel = 1.0 / del
	x2   = 1.0 - xo * rdel
	y2   = 1.0 - yo * rdel

c   'loca' contains the address where a block of data will be written.
c   a linked list is formed by 'next'.
	do 11 i = 1, nblk
	  ioff(i) = 1
11	  loca(i) = i
	next = nblk + 1
	do 12 i = 1, npb3 * nblk
12	  wrk(i) = dval

c  read data, find pigeon hole

100	go to ( 101, 102, 103, 104 ) itype
101	read( iswt, end=50,err=200,iostat=icheck ) x, y, z
	go to 110
102	read( iswt, end=50,err=200,iostat=icheck ) id, za
	x = za(ixchan)
	y = za(iychan)
	z = za(izchan)
	go to 110
103	read( iswt, ifmt, end=50,err=200,iostat=icheck ) x, y, z
	go to 110
104	read( iswt,    *, end=50,err=200,iostat=icheck ) x, y, z

110	if ( idirx .ne. 0 ) x = -x
	if ( x .lt. xmin  .or.  x .gt. xmax ) go to 100
	if ( y .lt. ymin  .or.  y .gt. ymax ) go to 100
	if ( z .lt. gmin  .or.  z .gt. gmax ) go to 100

	ib   = 0
	ic   = ic + 1
c     x&y converted to grid units
	x    = x * rdel + x2
	y    = y * rdel + y2
	bx   = x * dux
	by   = y * duy
	ibx  = int( bx + xlmt )
	iby  = int( by + ylmt )
	tstx = bx - float( ibx - 1 )
	tsty = by - float( iby - 1 )
	if ( ibx .gt. nsec ) then
	  ibx  = nsec
	  tstx = 1.0
	endif
	if ( iby .gt. ntier ) then
	  iby   = ntier
	  tsty  = 1.0
	endif
	ibset = nsec * ( iby - 1 ) + ibx
	mblk  = ibset

c  put data in pigeon hole, output when full
19	ixs  = npb3 * ( mblk - 1 )
	ip   = ixs  + ioff(mblk)
	wrk(ip)    = x
	wrk(ip+1)  = y
	wrk(ip+2)  = z
	ioff(mblk) = ioff(mblk) + 3
	if ( ioff(mblk) .ge. npb3 ) then
	  ndp        = ( ioff(mblk) - 1 ) / 3
	  call wrblk2( jswt, loca(mblk), next, ndp, nxyz, wrk(ixs+1) )
	  loca(mblk) = next
	  next       = next + 1
	  ioff(mblk) = 1
	endif
c
c  is data in overlap area ?
22	ib=ib+1
	go to (23,24,25,100)ib
23	if(tsty.gt.fy .or. iby.eq.1) go to 22
	mblk=ibset-nsec
	go to 19
24	if(ibx.eq.1 .or. iby.eq.1) go to 22
	if(tsty.gt.fy .or. tstx.gt.fx) go to 22
	mblk=ibset-nsec-1
	go to 19
25	if(tstx.gt.fx .or. ibx.eq.1) go to 100
	mblk=ibset-1
	go to 19

c   output unfilled pigeon holes
50	maxrec=next-1
	ip=1
	next=0
	do 51 i=1,nblk
	ndp=(ioff(i)-1)/3
	call wrblk2( jswt, loca(i), next, ndp, nxyz, wrk(ip) )
51	ip=ip+npb3

	print *, ic, '  data points in area'
	if ( ic .gt. 0 ) return
	print *,' data and grid coordinate mismatch ?'
       ierr = 1
       return
200    print *,' error no. ',icheck,' reading input data'
       ierr = 1
       return
	end
c******************************************************************************
	subroutine expdat( iswt, jswt, del, xo, yo, npb3,
     1                  buff, nword, ndmax, ntotal, ierr )

c  Input express file. 

	common /fmt/    ifmt(14), izchan, idirx, maxrec
	common /gparm/  mxc, mxr, nc, nr, nsec, ntier, lap, nim, epsm,
     1               dval, slope, gmin, gmax
	common /contn1/ loca(375), ioff(375)
	common /array/  wrk(40000)

	dimension       buff(nword,ndmax)

       ierr = 0

	ix = 1
	iy = 2
	iz = 3
	if ( izchan .ne. 0 ) then
c  default word assignments
c              id   channel
	  ix  = 2 +   1
	  iy  = 2 +   2
	  iz  = 2 +   izchan
	endif

	nxyz = npb3
	jrec = 0
	ic   = 0
	nblk = nsec * ntier
	nchk = npb3 - 3
	npb  = npb3 / 3
	err  = 1.e-2
	fudg = del * err
	xmax = float( mxc - 1 ) * del + xo - fudg
	ymax = float( mxr - 1 ) * del + yo - fudg
	xmin = xo + fudg
	ymin = yo + fudg

c  dimension of pigeon hole minus overlap
	cl   = float( nc - lap )
	rl   = float( nr - lap )
	dux  = 1. / cl
	duy  = 1. / rl
	fx   = ( float( lap ) - 0.1 ) * dux
	fy   = ( float( lap ) - 0.1 ) * duy

c  ensures smallest index (igx,y) in 'bwts' is 1
	err2 = 1.0 + 0.5 * err
	xlmt = 1.0 - err2 * dux
	ylmt = 1.0 - err2 * duy
	rdel = 1.0 / del
	x2   = 1.0 - xo * rdel
	y2   = 1.0 - yo * rdel

c   'loca' contains the address where a block of data will be written.
c   a linked list is formed by 'next'.
	do 10 i = 1, nblk
	  ioff(i) = 1
	  loca(i) = i
   10   continue
	do 20 i = 1, npb3 * nblk
	  wrk(i) = dval
   20   continue
	next = nblk + 1

c  fill buffer
50	idata = 0
	read( iswt, end=200 ) nw, ndata, nt, buff
	if ( nw .ne. nword  .or.  ndata .gt. ndmax  .or. 
     1    nt .ne. ntotal ) then
	  print *, ' nword, ndata, ntotal',   nw, ndata, nt
	  print *, ' nword, ndmax, ntotal',nword, ndmax, ntotal
	  print *, ' bad data file'
         ierr = 1
	  return
	endif

100	idata = idata + 1
	if ( idata .gt. ndata ) go to 50
	x = buff(ix,idata)
	y = buff(iy,idata)
	z = buff(iz,idata)

	if ( idirx .ne. 0 ) x = -x
	if ( x .lt. xmin  .or.  x .gt. xmax ) go to 100
	if ( y .lt. ymin  .or.  y .gt. ymax ) go to 100
	if ( z .lt. gmin  .or.  z .gt. gmax ) go to 100

	ib   = 0
	ic   = ic + 1
c     x&y converted to grid units
	x    = x * rdel + x2
	y    = y * rdel + y2
	bx   = x * dux
	by   = y * duy
	ibx  = int( bx + xlmt )
	iby  = int( by + ylmt )
	tstx = bx - float( ibx - 1 )
	tsty = by - float( iby - 1 )
	if ( ibx .gt. nsec ) then
	  ibx  = nsec
	  tstx = 1.0
	endif
	if ( iby .gt. ntier ) then
	  iby   = ntier
	  tsty  = 1.0
	endif
	ibset = nsec * ( iby - 1 ) + ibx
	mblk  = ibset

c  put data in pigeon hole, output when full
19	ixs  = npb3 * ( mblk - 1 )
	ip   = ixs  + ioff(mblk)
	wrk(ip)    = x
	wrk(ip+1)  = y
	wrk(ip+2)  = z
	ioff(mblk) = ioff(mblk) + 3
	if ( ioff(mblk) .ge. npb3 ) then
	  ndp        = ( ioff(mblk) - 1 ) / 3
	  call wrblk2( jswt, loca(mblk), next, ndp, nxyz, wrk(ixs+1) )
	  loca(mblk) = next
	  next       = next + 1
	  ioff(mblk) = 1
	endif

c  is data in overlap area ?
22	ib=ib+1
	go to (23,24,25,100)ib
23	if(tsty.gt.fy .or. iby.eq.1) go to 22
	mblk=ibset-nsec
	go to 19
24	if(ibx.eq.1 .or. iby.eq.1) go to 22
	if(tsty.gt.fy .or. tstx.gt.fx) go to 22
	mblk=ibset-nsec-1
	go to 19
25	if(tstx.gt.fx .or. ibx.eq.1) go to 100
	mblk=ibset-1
	go to 19
c
c   output unfilled pigeon holes
200	maxrec=next-1
	ip=1
	next=0
	do 51 i=1,nblk
	ndp=(ioff(i)-1)/3
	call wrblk2( jswt, loca(i), next, ndp, nxyz, wrk(ip) )
51	ip=ip+npb3
c
	print *, ic,'  data points in area'
	if ( ic .gt. 0 ) return
	print *,' data and grid coordinate mismatch ?'
       ierr = 1 
	return
	end
c******************************************************************************
	subroutine rejonl( region, xo, yo, del, xyz, npb3, ierr )

c  iterate regional grid coarse enough to fit into an array of length mxnode.

	parameter       ( mxnode=5000 )
	common /iodev/  igrid, iregn, imask, idata, irand, isubg
	common /contn1/ tmp(375), tmp2(375)
	common /fmt/    ifill(16), maxrec
	common /gparm/  mxc,mxr, nc,nr, nsec,ntier, lap, nim, epsm,
     1               dval ,slope1, gmin, gmax
	common /rgnl/   mc, mr, irxo, iryo, mdel, rg(mxnode)
	common /array/  zg(mxnode), iqd(mxnode), b(6*mxnode)

	dimension   irg(mxnode), xyz(npb3)
	equivalence ( irg, rg )
	character   id*56, p*8
	logical     lastt, lasts
	data nimr/30/, id/'regional grid'/, p/'minc'/, nz/1/
       ierr = 0
	npb = npb3 / 3

	mdel = 4
	mgu  = ( nsec  - 1 ) * ( nc - lap ) + nc - 1
	ngu  = ( ntier - 1 ) * ( nr - lap ) + nr - 1
10	mc   = mgu / mdel + 7
	mr   = ngu / mdel + 7
	nn   = mc * mr
	if ( nn .gt. mxnode ) then
	  mdel = mdel + 1
	  go to 10
	endif

	binit = 0.0
	do 15 i = 1, 6 * nn
15	  b(i) = binit

	do 16 i = 1, nn
	  rg(i)  = dval
16	  iqd(i) = 0

	slope = 1.0 / 2.0
	lastt = .false.
	dx    = nc - lap
	dy    = nr - lap
	endy  = dy
	rdel  = 1.0 / float( mdel )

	do 50 j=1,ntier
	endx=dx
	lasts=.false.
	if(j.eq.ntier) lastt=.true.
	do 40 i=1,nsec
	if(i.eq.nsec) lasts=.true.
	iadr=(j-1)*nsec+i
20	i2=1
	  if(iadr.gt.maxrec .or. iadr.le.0) then
	  print *,' rejonl: attempt to read record',iadr
         ierr = 1
	  return
	  endif
	read( irand, rec=iadr ) next, np, xyz
	iadr=next
	if(np.eq.0) go to 39
	if(np.gt.npb .or. np.lt.0) then
	  print *,'rejonl: np read error ',np
         ierr = 1
	  return
	  endif
	do 25 k=1,np*3-2,3
	k1=k+1
c   eliminate overlap in data
	if(lasts) go to 22
	if(xyz(k).gt.endx) go to 25
22	if(lastt) go to 23
	if(xyz(k1).gt.endy) go to 25
c   convert data to regional grid units
c   note that adding 3 puts the origin at -2*mdel
23	xyz(i2)=(xyz(k)-1.0)*rdel+3.0
	xyz(i2+1)=(xyz(k1)-1.0)*rdel+3.0
	xyz(i2+2)=xyz(k+2)
	i2=i2+3
25	continue
	np3=i2-1
	if(np3.lt.3) go to 39
	call assign(slope,0.,0.,mc,mr,b,np3,xyz)
39	if(next.ne.0) go to 20
40	endx=endx+dx
50	endy=endy+dy

	izset=1
	call bwts(izset,slope,binit,mc,mr,rg,iqd,b)
	call gridr(mc,mr,rg,tmp,dval,ierr)
       if(ierr.ne.0) return
	call curvmn(rg,iqd,b,mc,mr,epsm,nimr,st,end,ni)
	call glimit( mc, mr, rg, gmin, gmax, dval )

	irxo=1-2*mdel
	iryo=1-2*mdel
	if ( region .eq. 0.0 ) return

c  else write regional grid.
	dr=mdel*del
	rxo=xo-2.0*dr
	ryo=yo-2.0*dr
	write(iregn) id,p,mc,mr,nz,rxo,dr,ryo,dr
	k=1
	do 200 j=1,mr
	call rowio(mc,irg(k),0,iregn,iregn,ie)
200	k=k+mc
	return
	end
c******************************************************************************
	subroutine pcontl( nco, nro, ng, zg[huge], npb3, xyz, ierr )

	parameter       ( nwork=40000 )
	common /iodev/  igrid, iregn, imask, idata, irand, isubg
	common /gparm/  mxc, mxr, nc, nr, nsec, ntier, lap,
     1               nim, epsm, dval, slope, gmin, gmax
	common /qparm/  ihwind, mxcq, mxrq, maxq, maxg
	common /contn1/ zl(375), zb(375)
	common /assem/  ntot(2)
	common /array/  w(nwork)

	dimension iw(nwork), zg(ng)
	equivalence ( w(1), iw(1) )

	data nsav/3/, lout/6/

       ierr = 0

	ntot(1) = 0
	ntot(2) = 0
	maxg = ( nwork - ng ) / nco
	nn   = nwork / 8
	n2   = 1  + nn
	n3   = n2 + nn
	ir   = 0
	nbot = nc * nsav
	dx   = float( nc - lap )
	dy   = float( nr - lap )
	byo  = 0.0

	if ( ihwind .gt. 0 ) then
	  maxq = ( nwork - ng ) / mxcq
	  call padrow( ihwind, ntot(2), imask, mxcq, iw )
	endif

c  for each block, iterate a subgrid
	do 100 j = 1, ntier
	bxo = 0.0
	write( lout, 4 ) j
4	format(/,'   row section ',i2,/)

	do 50 i = 1, nsec
	if ( j .gt. 1 ) then
c  get lower boundary condition
c  caution zg(1) is w(1)
	  read( isubg, rec=i ) zg
	  i2 = nc * ( nr - lap ) + 1
	  do 5 ii = 1, nbot
	    zb(ii) = zg(i2)
	    i2     = i2 + 1
5	  continue
	endif

	mblk = nsec * ( j - 1 ) + i
	call icontl( i, j, ihwind, bxo, byo, xyz, npb3, ierr )
       if(ierr.ne.0) return
	if ( i .lt. nsec ) then
c  save leftside boundary
	  is = nc - lap + 1
	  i3 = 1
	  do 14 jj = 1, nr
	    i2 = is
	    do 12 ii = 1, nsav
	      zl(i3) = w(i2)
	      i3     = i3 + 1
	      i2     = i2 + 1
12	    continue
	    is = is + nc
14	  continue
	endif
	bxo = bxo + dx
50	continue

	call assemb( nc,   nr, iw(1),  nco, maxg, iw(ng+1),
     1             j,      0,  isubg, igrid, ierr )
       if(ierr.ne.0) return

	if ( ihwind .gt. 0 ) call assemb(
     1            nc,   nr, iw(1), mxcq, maxq, iw(ng+1),
     1             j, ihwind,  isubg, imask, ierr )
       if(ierr.ne.0) return

100	byo = byo + dy

	if ( ntot(1) .ne. nro ) then
	  print *, ' %% pcontl error: output grid contains'
	  print *, ' %%', ntot(1), 'rows instead of', nro
         ierr = 1
         return
	endif

	if ( ihwind .gt. 0 ) then
	  call padrow( ihwind, ntot(2), imask, mxcq, iw )
	  need = nro + 2 * ihwind
	  if ( ntot(2) .ne. need ) then
	    print *, ' %% pcontl error: masking grid contains'
	    print *, ' %%', ntot(2), 'rows instead of', need
           ierr = 1
           return
	  endif
	endif
	return
	end
c*****************************************
	subroutine padrow( nrow, jrow, iswt, nc, m[huge] )
	dimension m(nc)

	do 10 i = 1, nc
	  m(i) = 0
   10   continue
	do 20 i = 1, nrow
	  jrow = jrow + 1
	  call rowda( nc, m, 0, jrow, iswt, ierr )
   20   continue

	return
	end
c******************************************************************************
	subroutine icontl( ns, nt, ihwind, bxo, byo, xyz, npb3, ierr )

c  iteration control for one block
c  input xyz data, output iterated block

	parameter       ( mxnode=5000 )
	common /iodev/  igrid, iregn, imask, idata, irand, isubg
	common /gparm/  mxc,mxr, nc,nr, nsec,ntier, lap, nim, epsm,
     1               dval, slope, gmin, gmax
	common /contn1/ zl(375), zb(375)
	common /array/  zg(mxnode), iqd(mxnode), b(6*mxnode)

	dimension       izg(mxnode), xyz(npb3)
	equivalence     ( zg(1), izg(1) )
	data            lout/6/ 

       ierr = 0
	nn = nc * nr
	do 5 i=1,nn
5	iqd(i)=0
	call zgint(bxo,byo,ierr)
       if(ierr.ne.0) return
	binit=1.0
	if(slope.ne.0.0) binit=0.0
	do 10 i=1,nn*6
10	b(i)=binit

	idata = 0
	iadr  = nsec * ( nt - 1 ) + ns
20	read( irand, rec=iadr ) next, np, xyz
	np3   = 3 * np
	if ( np .ge. 1 ) then
	  idata = 1
	  call assign( slope, bxo, byo, nc, nr, b, np3, xyz )
	endif
	if ( next .eq. 0 ) go to 30
	iadr = next
	go to 20

30	if ( idata .gt. 0 ) then
	  izset = 0
	  call bwts( izset, slope, binit, nc, nr, zg, iqd, b )
	endif

	if ( ihwind .gt. 0 ) call wrblk( ns+nsec, iqd, nn, isubg )
	call contin( ns, nt, nc, nr )
	call gridr( nc, nr, zg, zl, dval, ier )
c  (zl used as work array)
	if ( ier .ne. 0 ) then
          print *, ' gridr error'
          ierr = 1
          return
       endif

c  apply minimum curvature equations.
	call curvmn( zg, iqd, b, nc, nr, epsm, nim, st, end, ni )

c  ad hoc limit to grid excursion.
	call glimit( nc, nr, zg, gmin, gmax, dval )

	write(lout,60) st,end,ni
60	format(' initial error=',1pe9.2,' end error=',e9.2,
     1      ' iterations=',i3)
       call wrblk( ns, izg, nn, isubg )
	return
	end
c******************************************************************************
	subroutine assign(slope,xogu,yogu,nc,nr,b[huge],np3,xyz)
c   associate np data values with grid locations of one block.
c   x,y,xogu,yogu coordinates are in grid units.
c   'slope' > 0 indicates that all data values within +-one half grid unit
c   will be combined by distance weighting.
c   'slope' = 0 indicates that the closest data value is the only one used.
	dimension xyz(np3),b(1)
	logical lslope
c
	if(np3.lt.3) return
	lslope=.false.
	if(slope.gt.0.0) lslope=.true.
c
	do 100 ixy=1,np3-2,3
	x=xyz(ixy)-xogu
	igx=int(x+.5)
	y=xyz(ixy+1)-yogu
	igy=int(y+.5)
c
	if(igx.gt.nc .or. igx.lt.1) go to 100
	if(igy.gt.nr .or. igy.lt.1) go to 100
c
30	z=xyz(ixy+2)
	ig=(igy-1)*nc+igx
	ib1=(ig-1)*6+1
	dx=x-float(igx)
	dy=y-float(igy)
	rsq=dx*dx+dy*dy
c
	if(lslope) go to 50
	if(rsq.ge.b(ib1)) go to 100
	b(ib1)=rsq
	b(ib1+1)=dx
	b(ib1+2)=dy
	b(ib1+5)=z
	go to 100
c
50	dwt=1.0/(rsq+slope)
	ib2=ib1+1
	ib3=ib1+2
	ib6=ib1+5
	b(ib1)=b(ib1)+dwt
	b(ib2)=b(ib2)+dx*dwt
	b(ib3)=b(ib3)+dy*dwt
	b(ib6)=b(ib6)+z*dwt
100	continue
	return
	end
c*****************************************
	subroutine glimit( nc, nr, zg, gmin, gmax, dval )
c  apply min/max limits to grid array.
	dimension zg(nc*nr)
	logical tstmin, tstmax

	tstmin = .false.
	tstmax = .false.
	if ( gmin .ne. -dval ) tstmin = .true.
	if ( gmax .ne.  dval ) tstmax = .true.
	if ( .not. ( tstmin .or. tstmax )  ) return

	do 10 i = 1, nc * nr
	  if ( tstmin  .and.  zg(i) .lt. gmin )  then
	    zg(i) = gmin
	    go to 10
	  endif
	  if ( tstmax  .and.  zg(i) .gt. gmax )  zg(i) = gmax
10	continue

	return
	end
c******************************************************************************
	subroutine bwts(izset,slope,binit,nc,nr,zg,iqd,b[huge])
c   initialize iqd and zg arrays and
c   calculate weights for minimum curvature equations.
c   'b' is dimensioned (6,nc*nr)
c   b(1,ig) contains either r squared or 1./weight
	dimension itabl(4)
	dimension zg(1),iqd(1),b(1)
	logical notslp
	data itabl/3,4,2,1/
c
	notslp=.true.
	if(slope.gt.0.0) notslp=.false.
	ib1=1
	do 100 ig=1,nc*nr
	if(b(ib1).eq.binit) go to 100
	ib2=ib1+1
	ib3=ib1+2
	ib6=ib1+5
	if(notslp) go to 50
	rwt=1./(b(ib1))
	b(ib2)=b(ib2)*rwt
	b(ib3)=b(ib3)*rwt
	b(ib6)=b(ib6)*rwt
	b(ib1)=b(ib2)*b(ib2)+b(ib3)*b(ib3)
c
50	if(b(ib1).gt.0.0025) go to 80
	iqd(ig)=-1
	zg(ig)=b(ib6)
	go to 100
80	ix=1
	if(b(ib2).ge.0.0) ix=2
	iy=0
	if(b(ib3).ge.0.0) iy=2
	iqd(ig)=itabl(ix+iy)
	if(izset.ne.0) zg(ig)=b(ib6)
100	ib1=ib1+6
c
c  calculate weighting
200	is=2*nc+3
	ie=3*nc-2
	do 230 j=3,nr-2
	ib1=(is-1)*6+1
	do 220 i=is,ie
	if(iqd(i))220,220,210
210	ib2=ib1+1
	ib3=ib2+1
	ib4=ib3+1
	ib5=ib4+1
	ib6=ib5+1
	dx=abs(b(ib2))
	dy=abs(b(ib3))
	dy2=dy*dy
	f1=dx*(dx+dy+dy+1.)
	b5=4./(f1+dy2+dy)
	b4=(b5*f1*.5)-1.
	b5dx=b5*dx
	b4b4=b4+b4
	b3=b5dx*(dy+1.)-b4b4
	b(ib2)=2.+b5dx-(b5*dy2+b4b4+b3)
	b(ib1)=b3+b4-b5dx
	b(ib3)=b3
	b(ib4)=b4
	b(ib5)=b5*b(ib6)
	b(ib6)=1./(1.+b(ib1)+b(ib2)+b3+b4+b5)
220	ib1=ib1+6
	is=is+nc
230	ie=ie+nc
	return
	end
c******************************************************************************
	subroutine zgint(bxo,byo,ierr)
	common /rgnl/ mc,mr,irxo,iryo,mdel,rg(5000)
	common /array/  zg(5000),iqd(5000),b(30000)
	common /gparm/  mxc,mxr,nc,nr,nsec,ntier,lap,
     1               nim, epsm, dv, slope, gmin, gmax
c
       ierr = 0
	igxs=1+int(bxo+.001)
	igys=1+int(byo+.001)
	igxe=nc+int(bxo+.001)
	igye=nr+int(byo+.001)
	irxs=1+(igxs-irxo)/mdel
	irys=1+(igys-iryo)/mdel
	irxe=1+(igxe-irxo)/mdel
	irye=1+(igye-iryo)/mdel
	if(irxs.lt.2 .or. irxe.gt.mc-2) then
          print *,' error 444'
          ierr = 1
          return
       endif
	if(irys.lt.2 .or. irye.gt.mr-2) then 
          print *,' error 555'
          ierr = 1
          return
       endif
	ncb = 1 + mdel*(irxe-irxs+1)
	nrb = 1 + mdel*(irye-irys+1)
	if(ncb*nrb.gt.30000) then
         print *, ' insert: ncb*nrb > 30,000'
         ierr = 1
         return
       endif
c
	do 100 i=1,ncb*nrb
100	b(i)=dv
	jb=1
	do 120 jr=irys,irye
	ib=1
	do 110 ir=irxs,irxe
	nint=mdel
	if(ir.eq.irxe .or. jr.eq.irye) nint=mdel+1
	call cell3(nint,mdel,ir,jr,mc,rg,zg)
	call insert(ib,jb,zg,nint,nint,b,ncb,nrb,ierr)
       if(ierr.ne.0) return
110	ib=ib+mdel
120	jb=jb+mdel
c
	do 200 i=1,nc*nr
200	zg(i)=dv
	ibs=irxo+mdel*(irxs-1)
	jbs=iryo+mdel*(irys-1)
	izgs=igxs-ibs+1
	jzgs=igys-jbs+1
	if(izgs.le.0 .or. jzgs.le.0) then 
          print *,' error 666'
          ierr = 1
          return
       endif
	izge=izgs-1+nc
	jzge=jzgs-1+nr
	if(izge.gt.ncb .or. jzge.gt.nrb) then 
          print *,' error 777'
          ierr = 1
          return
       endif
	do 210 jg=1,nr
	ij1=izgs+ncb*(jzgs+jg-2)
	ij2=1+nc*(jg-1)
	do 210 ig=1,nc
	zg(ij2)=b(ij1)
	ij1=ij1+1
210	ij2=ij2+1
c
	idv=0
	do 250 i=1,nc*nr
250	if(zg(i).eq.dv) idv=1
	if(idv.eq.1)then
         print *,' zgint: dvals present in zg'
         ierr = 1
       endif
	return
	end
c******************************************************************************
	subroutine cell3(nint,nsub,igx,igy,nc,zg,zint)
c  cubic interpolation first in y then in x.  
c  returns the interpolated z values for the interior of a cell
c  with lower-left address (igx,igy).
	dimension z(3),a1(4),a2(4),a3(4),zy(4)
	dimension zint(nint,nint),zg(1)
	del =1.0/float(nsub)
	ij=igx+nc*(igy-1)
c
	igs=ij-(nc+1)
	do 20 icol=1,4
	ig=igs
	z0=zg(ig)
	do 10 i=1,3
	ig=ig+nc
10	z(i)=zg(ig)-z0
	a1(icol)=  3.0*z(1) - 1.5*z(2) + .3333333*z(3)
	a2(icol)= -2.5*z(1) + 2.0*z(2) -      0.5*z(3)
	a3(icol)=  0.5*z(1) - 0.5*z(2) + .1666667*z(3)
20	igs=igs+1
c
	y=1.0
	do 100 jout=1,nint
	ig=ij-(nc+1)
	do 30 icol=1,4
	zy(icol)= zg(ig) + y*(a1(icol) + y*(a2(icol) + y*a3(icol)))
30	ig=ig+1
	do 40 i=1,3
40	z(i)=zy(i+1)-zy(1)
	b1 =  3.0*z(1) - 1.5*z(2) + .3333333*z(3)
	b2 = -2.5*z(1) + 2.0*z(2) -      0.5*z(3)
	b3 =  0.5*z(1) - 0.5*z(2) + .1666667*z(3)
	x=1.0
	do 50 iout=1,nint
	zint(iout,jout)= zy(1) + x*(b1 + x*(b2 + x*b3))
50	x=x+del
100	y=y+del
	return
	end
c******************************************************************************
	subroutine insert(ib,jb,a,ma,na,b[huge],mb,nb,ierr)
c  insert embedded matrix 'a' into 'b'
	dimension a(ma,na), b(mb,nb)
	ierr=0
	if(ib+ma-1.gt.mb) ierr=1
	if(jb+na-1.gt.nb) ierr=ierr+10
	if(ierr.eq.0) go to 5
	print *,' insert: ierr=',ierr
	print *,ib,jb,ma,na,mb,nb
       ierr = 1
	return
5	j2=jb
	do 20 j=1,na
	i2=ib
	do 10 i=1,ma
	b(i2,j2)=a(i,j)
10	i2=i2+1
20	j2=j2+1
	return
	end
c******************************************************************************
	subroutine contin(ns,nt,nc,nr)
	common /array/  zg(5000),iqd(5000),b(30000)
	common /contn1/ zl(375),zb(375)
	nsav=3
	if(ns.eq.1) go to 20
c  insert leftside boundary condition
	j=1
	igs=1
	do 15 ii=1,nr
	ig=igs
	do 10 i=1,nsav
	zg(ig)=zl(j)
	if(i.ne.nsav) iqd(ig)=-1
	ig=ig+1
10	j=j+1
15	igs=igs+nc
c
20	if(nt.eq.1) go to 30
c  insert bottom boundary condition
	do 24 i=1,2*nc
	zg(i)=zb(i)
24	iqd(i)=-1
	do 25 i=2*nc+1,3*nc
25	zg(i)=zb(i)
30	return
	end
c******************************************************************************
	subroutine gridr(nc,nr,zg,wz,dval,ierr)
c  initialize grid with reasonable anomalies
c  dval is the flag value to be replaced
c  wz at least max(nc,nr)
	dimension iw1(3)
	dimension zg(1),wz(1)
	data iw1/3,7,9/
	ierr=0
	nn=nc*nr
	nsep=(iw1(3)-1)/2
	ns1=nsep+1
	mxpass=max(nc,nr)/(2*nsep)
	i=1
201	if(i.gt.nn) return
	if(zg(i).ge.dval) go to 200
	i=i+1
	go to 201
c
c  insert control point net
c     distances for ring averages specified by iw1
200	if(nc.lt.nsep+ns1 .or. nr.lt.nsep+ns1) go to 140
	ipass=0
	iflag=1
130	iw=iw1(iflag)
	ihw=(iw-1)/2
	iset=ihw*nc+ihw
131	nass=0
	do 110 jj=ns1,nr-ihw,nsep
	ip=(jj-1)*nc+ns1
	do 100 ii=ns1,nc-ihw,nsep
	if(zg(ip).lt.dval) go to 100
	ips=ip-iset
	it=0
	t=0.0
	do 121 j=1,iw
	ip2=ips
	do 120 i=1,iw
	if(zg(ip2).ge.dval) go to 120
	t=zg(ip2)+t
	it=it+1
120	ip2=ip2+1
121	ips=ips+nc
	if(it-3.lt.0) go to 100
	zg(ip)=t/float(it)
	nass=nass+1
100	ip=ip+nsep
110	continue
	ipass=ipass+1
	if(nass.eq.0 .and. iflag.eq.3) go to 140
	if(ipass.gt.mxpass) go to 140
	if(iflag-2)133,132,131
133	iflag=2
	go to 130
132	iflag=3
	go to 130
140	continue
c
c  fill holes
	inc=nc*nsep
	j=inc+1
	if(nr.lt.ns1) go to 21
	do 20 irow=ns1,nr-nsep,nsep
	call plugm(nc,zg(j),dval)
20	j=j+inc
21	do 28 icol=1,nc
	j=icol
	do 22 k=1,nr
	wz(k)=zg(j)
22	j=j+nc
	call plugm(nr,wz,dval)
	j=icol
	do 24 k=1,nr
	zg(j)=wz(k)
24	j=j+nc
28	continue
c  final check
	do 40 i=1,nn
	if(zg(i).ge.dval) go to 41
40	continue
	return
41	t=0.0
	it=0
	do 42 i=1,nn
	if(zg(i).ge.dval) go to 42
	t=t+zg(i)
	it=it+1
42	continue
	if(it.eq.0) then
         print *, ' cannot init grid'
         ierr = 1
         return
       endif
	t=t/float(it)
	print  43,t
43	format(' gridr init with',1pe15.5)
	do 44 i=1,nn
	if(zg(i).ge.dval) zg(i)=t
44	continue
	return
	end
c******************************************************************************
	subroutine curvmn(zg,iqd,b[huge],nc,nr,epsmx,nim,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  Mike Webring, USGS Open-File report 81-1224.
	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 assemb( nc, nr, m[huge], ncout, nmax, mw[huge], 
     1                  nt, ihwind, iswt, jswt, ierr )

c  input a completed tier of blocks and assemble row records
c  for output as the finished grid.
c  ihwind > 0 is the switch indicating masking row is produced.

	common /gparm/  mxc, mxr, idum(2), nsec, ntier, lap,
     1               nim, epsm, dval, slope, gmin, gmax
	common /assem/  ntot(2)
	dimension       m(nc,nr), mw(ncout,nmax)

       ierr = 0
	nco  = mxc - 4
	nro  = mxr - 4
	nc1  = nc  - lap
	nr1  = nr  - lap

	ioff = 0
	if ( ihwind .gt. 0 ) ioff = nsec

	if ( nt .eq. 1 ) then
	  js  = 3
	  nt1 = 999
	  else
	  js  = 1
	  nt1 = nt
	endif

	nrout = nro
	if ( ntier .gt. 1 ) then
	  ntest = ntier - nt1
	  nrout = nr1   - 2
	  if ( ntest .gt. 0 ) nrout = nr1
	  if ( ntest .eq. 0 ) nrout = nro - nr1 * ( ntier - 1 ) + 2
	endif
	if ( nrout .gt. nr ) then
          print *, ' %%assemb: nrout error'
          ierr = 1
          return
       endif
	ntime = 0
4	ntime = ntime + 1
	if ( nrout - ntime * nmax  .gt.  0 ) go to 4

	nrow = nrout
	if ( nrout .gt. nmax ) nrow = nrout / ntime + 1

	do 200 itime = 1, ntime
	istop = nco + 2
	itot  = 0
	is    = 3
	ie    = nc1
	iws   = ihwind + 1
	if ( itime .eq. ntime ) nrow = nrout - nrow * ( itime - 1 )

	do 10 isec = 1, nsec
	  irec = isec + ioff
	  read( iswt, rec=irec ) m

	  if ( ihwind .gt. 0 ) then
c  set data flags in masking grid.
	    do 5 j = 1, nr
	      do 5 i = 1, nc
	      if ( m(i,j) .ne. 0 ) m(i,j) = 1
5	    continue
	  endif

	  jin = js
	  if ( isec .eq. nsec ) ie = istop

	  do 8 j = 1, nrow
	    iwcol = iws
	    do 7 i = is, ie
	      mw(iwcol,j) = m(i,jin)
7	      iwcol       = iwcol + 1
	    jin=jin+1
8	  continue

	  itot  = ( ie - is + 1 ) + itot
	  istop = nco - itot
	  iws   = iwcol
	  is    = 1
10	continue

	do 100 j = 1, nrow
	  if ( ihwind .gt. 0 ) then
c  write masking row
	    i2 = ncout
	    do 50 i = 1, ihwind
	      mw(i,j)  = 0
	      mw(i2,j) = 0
	      i2       = i2 - 1
50	    continue
	    ntot(2) = ntot(2) + 1
           call rowda( ncout, mw(1,j), 0, ntot(2), jswt, ierr )
	    if ( ntot(2) .eq. nro + ihwind ) go to 999
	    else
c  write data row
	    ntot(1) = ntot(1) + 1
           call rowio( ncout, mw(1,j), 0, jswt, jswt, ierr )
	    if ( ntot(1) .eq. nro ) go to 999
	  endif
100	continue

	  js = js + nrow
200	continue

999	return
	end
c******************************************************************************
	subroutine radmsk( np, iswt, ncol, iwind,
     &                    mbuff[huge], msum[huge], mout[huge],ierr )
c  generate masking grid one radius around all nozero input grid points.
c  input random access grid and a radius.
c  output masking grid of 0 and 1's.

	dimension   mbuff(ncol,iwind), msum(ncol), mout(ncol)
	character   title*56, pgm*8

       ierr =0
	ihwind = iwind / 2
	iwindw = 2 * ihwind + 1
	if ( iwindw .gt. iwind ) then
	  ihwind = ihwind - 1
	  iwindw = 2 * ihwind + 1
	endif
	npmin  = np
	if ( npmin .lt. 1 ) npmin = 1
	if ( npmin .gt. iwindw * iwindw ) npmin = iwindw * iwindw

	read( iswt, rec=1 ) title, pgm, nc, nr
	if ( ncol .lt. nc ) then
	  print *, ' %%radmsk error: not enough working area'
         ierr = 1
	  go to 999
	endif
	if ( ihwind .lt. 1 ) then
	  print *, ' %%radmsk error: window less than 3 cells'
         ierr = 1
	  go to 999
	endif
	istop = nc - ihwind
	jstop = nr - ihwind

c  initialize buffer
	do 1 j = 1, iwindw
	  irow = j
	  call rowmsk( nc, mbuff(1,j), irow, iswt, ierr )
1	continue
	do 2 i = 1, nc
	  msum(i) = 0
2	continue
	do 3 j = 1, iwindw
	  do 3 i = 1, nc
	    msum(i) = msum(i) + mbuff(i,j)
3	continue

	jptr = 1
	irow = iwindw
	do 100 jrow = ihwind + 1, jstop

	call rowmsk( nc, mout, jrow, iswt, ierr )

	nwind = 0
	do 20 i = 1, iwindw - 1
20	  nwind = nwind + msum(i)
	iptr1 = 1
	iptr2 = iwindw

	do 30 i = ihwind + 1, istop
	  nwind = nwind + msum(iptr2)
	  if ( nwind .ge. npmin ) mout(i) = 1
	  nwind = nwind - msum(iptr1)
	  iptr1 = iptr1 + 1
	  iptr2 = iptr2 + 1
30	continue

	call rowda( nc, mout, 0, jrow, iswt, ierr )
	if ( jrow .eq. jstop ) go to 999

	do 40 i = 1, nc
40	  msum(i) = msum(i) - mbuff(i,jptr)
	irow = irow + 1
	if ( irow .gt. nr ) go to 999
	call rowmsk( nc, mbuff(1,jptr), irow, iswt, ierr )
	if ( ierr .ne. 0 ) go to 999
	do 50 i = 1, nc
50	  msum(i) = msum(i) + mbuff(i,jptr)

	jptr = jptr + 1
	if ( jptr .gt. iwindw ) jptr = 1
100	continue

999	if ( jrow .ne. nr - ihwind ) then
	  print *, ' %%radmsk warning: ihwind, nr, irow, jrow='
	  print *,                     ihwind, nr, irow, jrow
         ierr = 1
	endif
	return
	end
c*******************************************
	subroutine outside( imask, ihwind, iwork[huge], nwork, ierr )

c  fill the polygons of no-data that intersect the grid boundary.

	dimension iwork(nwork)
	character title*56, pgm*8

	read( imask, rec=1 ) title, pgm, ncm, nrm
	call setpoly( imask, ihwind, iwork, ncm, nrm, ierr )
       if(ierr.ne.0) return

	nrblk = nwork / ncm
	call fill( imask, iwork, ncm, nrm, nrblk )

	return
	end
c*******************************************
	subroutine setpoly( iswt, ihwind, m[huge], nc, nr,ierr )
	dimension m(nc)
	character title*56, pgm*8

c  set no-data sections on boundary to seed of 2.
c  flag half window section with 3.

c  col&row locations for the boundary of the output grid.
       ierr = 0
	icolm = ihwind + 1
	irowm = ihwind + 1
	icolx = nc - ihwind
	irowx = nr - ihwind

c  bottom half window
	do 20 j = 1, irowm - 1
	  irec = j + 1
	  read( iswt, rec=irec ) iy, m
	  do 10 i = 1, nc
	    m(i) = 3
   10     continue
	  write( iswt, rec=irec ) iy, m
   20   continue

c  middle section 

	do 60 j = irowm, irowx
	  irec = j + 1
	  read( iswt, rec=irec ) iy, m
c  left&right half window
	  i2 = nc
	  do 40 i = 1, ihwind
	    m(i)  = 3
	    m(i2) = 3
	    i2    = i2 - 1
   40     continue
c  enter seed values on output grid boundary.
	  if ( m(icolm) .eq. 0 ) m(icolm) = 2
	  if ( m(icolx) .eq. 0 ) m(icolx) = 2
	  if ( j .eq. irowm  .or.  j .eq. irowx ) then
	    do 50 i = icolm, icolx
	      if ( m(i) .eq. 0 ) m(i) = 2
   50       continue
	  endif
	  write( iswt, rec=irec ) iy, m
   60   continue

c  top half window
	do 80 j = irowx + 1, nr
	  irec = j + 1
	  read( iswt, rec=irec, err=200 ) iy, m
	  do 70 i = 1, nc
	    m(i) = 3
   70     continue
	  write( iswt, rec=irec ) iy, m
   80   continue

	return
200	read( iswt, rec=1 ) title, pgm, ncg, nrg
	print *, ' %%setpoly: EOF at row', j, ' of', nrg
       ierr = 1
	return
	end
c*******************************************
		subroutine fill( iswt, m[huge], nc, nrmax, nrbuf )
c  seed fill driver. 
	dimension js(20), je(20), m(nc,nrbuf)

c  setup buffer dimensions
	jbuf  = 1
	jdir  = 1
	nbuf  = 1
	js(1) = 1
	je(1) = nrmax
	if ( nrmax .le. nrbuf ) go to 20

c  each buffer will share its ending row with the start of the next.
	nbuf = nrmax / nrbuf + 1
 5	nr   = ( nrmax + nbuf ) / nbuf + 1
	je(1) = nr
	do 10 i = 2, nbuf
	  js(i) = je(i-1)
	  je(i) = js(i) + nr - 1
	  if ( je(i) .ge. nrmax ) then
	    je(i) = nrmax
	    if ( i .ne. nbuf ) nbuf = i
	    go to 20
	  endif
   10   continue
c  increase number of buffers
	nbuf = nbuf + 1
	go to 5  

c  cycle buffers

20	nr   = je(jbuf) - js(jbuf) + 1
	jrow = js(jbuf)
	do 30 j = 1, nr
	  call rowda( nc, m(1,j), -1, jrow, iswt, ierr )
	  jrow = jrow + 1
   30   continue

	call fillblk( m, nc, nr, nfill ) 
	ndir = ndir + nfill

	if ( nfill .gt. 0 ) then
	  jrow = js(jbuf)
	  do 40 j = 1, nr
	    call rowda( nc, m(1,j), 0, jrow, iswt, ierr )
	    jrow = jrow + 1
   40     continue
	endif

	if ( nbuf .eq. 1 ) go to 999
	if ( ( jdir .gt. 0  .and.  jbuf .eq. nbuf )  .or. 
     1    ( jdir .lt. 0  .and.  jbuf .eq. 1 ) )   then
c  pass completed 
	  if ( ndir .eq. 0 ) then
	    if ( noass .eq. 1 ) go to 999
	    noass = 1
	    else
	    noass = 0
	  endif
	  jdir = -jdir
	  ndir = 0
	endif
	
	jbuf = jbuf + jdir
	go to 20

999	return
	end
c*******************************************
		subroutine fillblk( m[huge], nc, nr, nfill )
c  seed fill controller for one block.
c  at least one element of m should equal 2.
	common /filcom/ jcurr, jlast
	dimension       m(nc,nr)
	logical         nearby

	nfill = 0
	noass = 0
	npass = 0

	jmin  = 1
	jmax  = nr
	jmina =  9999
	jmaxa = -9999
	nearby = .false.

10	npass = npass + 1
	ndir  = 0

c  scan up
	jcurr = jmin
   15   if ( nearby  .or.  jcurr .le. jmax + 1 ) then
	  if ( jcurr .gt. nr ) go to 20
	  if ( jcurr .lt. 2  ) jcurr = 2
	  nass  = 0
	  jlast = jcurr - 1
	  call fillrow( nc, m(1,jcurr), m(1,jlast), nass )
	  if ( nass .gt. 0 ) then
	    nearby = .true. 
	    ndir   = ndir + nass
	    if ( jcurr .gt. jmaxa ) jmaxa = jcurr
	    else
	    nearby = .false.
	  endif
	  jcurr = jcurr + 1
        go to 15
        endif

c  scan down
20	jcurr = jmax
   25   if ( nearby  .or.  jcurr .ge. jmin - 1 ) then
	  if ( jcurr .lt. 1 ) go to 30
	  if ( jcurr .gt. nr - 1 ) jcurr = nr - 1
	  nass  = 0
	  jlast = jcurr + 1
	  call fillrow( nc, m(1,jcurr), m(1,jlast), nass )
	  if ( nass .gt. 0 ) then
	    nearby = .true. 
	    ndir   = ndir + nass
	    if ( jcurr .lt. jmina ) jmina = jcurr
	    else
	    nearby = .false.
	  endif
	  jcurr = jcurr - 1
        go to 25
        endif

30	nfill = nfill + ndir

	if ( ndir .eq. 0 ) then
	  if ( noass .eq. 1 ) go to 999
	  noass = 1
	  jmin  = 1
	  jmax  = nr
	  else
	  noass = 0
	  jmin  = jmina
	  jmax  = jmaxa
	endif
	go to 10

999	return
	end
c*******************************************
		subroutine fillrow ( nc, mcurr[huge], mlast[huge], nass )
c  scan row for seed values.
	common /filcom/ jcurr, jlast
	dimension mcurr(nc), mlast(nc)

	do 100 icurr = 2, nc - 1
	  if ( mcurr(icurr) .eq. 0 ) then

	    if ( mcurr(icurr-1) .eq. 2 ) then
	      mcurr(icurr) = 2
	      nass         = nass + 1
	      go to 100
	    endif

	    if ( mlast(icurr) .eq. 2 ) then
	      mcurr(icurr) = 2
	      nass         = nass + 1
	      if ( mcurr(icurr-1) .eq. 0 ) 
     1          call backfill( mcurr, icurr-1, nass )
	      go to 100
	    endif

	    if ( mcurr(icurr+1) .eq. 2 ) 
     1        call backfill( mcurr, icurr, nass )
	  endif
100	continue

	if ( mcurr(1) .eq. 0  .and. 
     1  ( mcurr(2) .eq. 2  .or.  mlast(1) .eq. 2 ) ) then
	  mcurr(1) = 2
	  nass     = nass + 1
	endif

	if ( mcurr(nc)   .eq. 0  .and.
     1  ( mcurr(nc-1) .eq. 2  .or.  mlast(nc) .eq. 2 ) ) then
	  mcurr(nc) = 2
	  nass      = nass + 1
	endif

	return
	end
c*******************************************
	subroutine backfill( mcurr[huge], icurr, nass ) 
	dimension mcurr(icurr)
	do 10 i = icurr, 1, -1
	  if( mcurr(i) .ne. 0 ) return
	  mcurr(i) = 2
	  nass     = nass + 1
   10   continue
	return
	end
c*******************************************
	subroutine rowmsk( nc, mask[huge],  irow, iswt, ierr )
	dimension mask(nc)
	ierr = 0
	irec = irow + 1
	read( iswt, rec=irec, err=20 ) iy, mask
	do 10 i = 1, nc
	  if ( mask(i) .ne. 0 ) mask(i) = 1
10	continue	    
	return
20	ierr = 1
	return
	end
c*******************************************
	subroutine templat( istart, jstart, iflag, dval,
     1                   iswt, m[huge], ncm, nrm,
     1                   jswt, z[huge], nc,  nr, ierr   )

c iswt contains the template where dval's are inserted into jswt.
c jswt area must be a subset of iswt.

	dimension m(ncm), z(nc)
       ierr = 0
	if ( istart .le. 0 ) go to 200
	if ( jstart .le. 0 ) go to 200
	ioff = istart - 1
	joff = jstart - 1
	iend = ioff + nc
	jend = joff + nr
	if ( iend .gt. ncm ) go to 200
	if ( jend .gt. nrm ) go to 200

	do 20 j = 1, nr
c  grid header is record 1, row 1 is record 2.
	  irec = ( j + 1 ) + joff
c        jrec = ( j + 1 )
	  read( iswt, rec=irec ) iy, m
	  read( jswt ) y,  z

	  im = 1 + ioff
	  do 10 i = 1, nc
	    if ( m(im) .eq. iflag ) z(i) = dval
	    im = im + 1
   10     continue

	  write( jswt+10 ) y, z
   20   continue

	return
200	print *, ' %%templat: grid is not a subset of the template.'
       ierr = 1
	return
	end
c******************************************************************************
	subroutine plugm(n,z,dv)
c  plug holes using linear interpolation
	dimension z(n)
	do 1 is=1,n
	if(z(is) .ne. dv) go to 2
1	continue
	return
2	ix=is
3	ix=ix-1
	if(ix.lt.1) go to 4
	z(ix)=z(is)
	go to 3
4	do 5 idv=is,n
	if(z(idv) .eq. dv) go to 6
5	continue
	return
6	is=idv-1
	do 7 ie=idv,n
	if(z(ie) .ne. dv) go to 10
7	continue
	ix=is
9	ix=ix+1
	if(ix.gt.n) return
	z(ix)=z(is)
	go to 9
10	dz=(z(ie)-z(is))/float(ie-is)
	do 11 i=is+1,ie-1
11	z(i)=z(i-1)+dz
	is=ie
	go to 4
	end
	subroutine rowio(n,iz[huge],iop,idev,jdev,iend)
c  where read iop<0; write iop=0; r&w iop>0
	dimension iz(n)
	y=0.
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,iz
	if(iop)9,9,2
2	write(jdev) y,iz
9	return
10	iend=1
	return
	end
c******************************************************************************
	subroutine rowda( n, iz[huge], iop, irow, idev, ierr )
c  header is record 1, so irec = irow + 1
	dimension iz(n)
	ierr = 0
	irec = irow + 1
	if ( iop .eq. 0 ) then
	  write( idev, rec=irec, iostat=ierr ) ydum, iz
	  else
	  read ( idev, rec=irec, iostat=ierr ) ydum, iz
	endif
	return
	end
c*******************************************
	subroutine wrblk( loc, m, n, iswt )
	dimension m(n)
	write( iswt, rec=loc ) m
	return
	end
c*******************************************
	subroutine wrblk2( iswt, irec, next, ndp, nxyz, xyz[huge] )
	dimension xyz(nxyz)
	write( iswt, rec=irec ) next, ndp, xyz
	return
	end
c******************************************************************************
      subroutine namemc(icmd,ierr)
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(numt=50)
      parameter(ivar=80)
      character*80 tinp(numt)
      character*6 tindex(numt)
      integer*2 ix
      logical array
      common/lout/tindex,ix,tinp,array
      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/' '/,sls/'/'/
c
c     get namelist start name(e.g.&parms)
c
      ierr=0
      ix=0
      it=0
    5 if(icmd.ne.-99) then
        read(icmd,1000,end=910) var
      else
        go to 10
      endif
      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 if(icmd.ne.-99) then
        read(icmd,1000,end=900) var
 1000 format(80a1)
      else
        it = it + 1
        read(tinp(it),'(a)') var
      endif
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar.or.var(i)
     & .eq.sls) 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
      call clearscreen($GCLEARSCREEN)
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      ierr=1
      return
      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
     & .or.var(l).eq.sls) 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
      call clearscreen($GCLEARSCREEN)
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      ierr=1
      return
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      call clearscreen($GCLEARSCREEN)
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      ierr=1
      return
      endif
   80 continue
c     no value on same line as variable
c     read next line
      if(icmd.ne.-99) then
        read(icmd,1000,end=900) var
      else
        it = it + 1
        read(tinp(it),'(a)') var
      endif
      j=0
      go to 50
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum,ierr)
      if(ierr.ne.0) return
      tvar=blank
      k=l+1
      if(k.gt.80) go to 10
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
     & .or.var(l).eq.sls).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.or.
     & var(l).eq.sls) 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
      if(nn.ne.0) then
        inum=inum+1
        go to 90
      endif
      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.or.
     & var(j).eq.sls) 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
      call clearscreen($GCLEARSCREEN)
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      ierr=1
      return
      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.or.
     & var(l).eq.sls) 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 $'
      ierr=1
      return
      end
c*******************************************************************
      subroutine check(pvar,tvar,nn,chv,nvar,numa,inum,ierr)
c
c     assigns values to proper variable
c     variables are passed to program mcgrd 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 mcgrd
c
      parameter(nnvar=21,numr=7)
      character*6 pvar,var(nnvar)
      character*56 tvar,kvar,cfmt
      logical chv
      character*56 title,ifile,ofile,ifmt
      common/fmt/ifmt,izchan,idirx,maxrec
      common/gparm/mxc,mxr,nc,nr,nsec,ntier,lapovr,
     & nim,epsm,dv,slope,gmin,gmax
      common/gparm2/cm(3),baslat(3)
      common/nlparm/ifile,ofile,title,xo,yo,del,nco,
     & nro,npmin,radius,hull,whole,region,iproj,cmdeg,bldeg
      data var/'nc','nr','idirx','nim','npmin','izchan',
     & 'xo','yo','del','epsm','whole','radius','region','slope',
     & 'ifile','ofile','title','ifmt','hull','gmin','gmax'/
      numa=22
      if(nn.gt.56) then
        print *,'more than 56 characters in variable: ',pvar
        ierr = 1
        return
      endif
      do 190 i=1,nnvar
      if(pvar.ne.var(i)) go to 190
      if(.not.chv) then
c
c     noncharacter value
c     right justify then 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
      go to (101,102,103,104,106,107,108,109,110,111,112,
     & 113,114,115,116,117,118,119,120,121,122),i
  101 nc=jvar
      go to 200
  102 nr=jvar
      go to 200
  103 idirx=jvar
      go to 200
  104 nim=jvar
      go to 200
  106 npmin=jvar
      go to 200
  107 izchan=jvar
      go to 200
  108 xo=xvar
      go to 200
  109 yo=xvar
      go to 200
  110 del=xvar
      go to 200
  111 epsm=xvar
      go to 200
  112 whole=xvar
      go to 200
  113 radius=xvar
      go to 200
  114 region=xvar
      go to 200
  115 slope=xvar
      go to 200
  116 ifile=tvar(1:nn)
      go to 200
  117 ofile=tvar(1:nn)
      go to 200
  118 title=tvar(1:nn)
      go to 200
  119 ifmt=tvar(1:nn)
      go to 200
  120 hull=xvar
      go to 200
  121 gmin=xvar
      go to 200
  122 gmax=xvar
      go to 200
  190 continue
      write(*,*) ' error in namelist - ',pvar,' variable not included'
      ierr=1
      return
  200 nvar=i
      call tcheck(pvar,tvar,nn,chv)
      return
      end
c***********************************************************************
      subroutine tcheck(pvar,tvar,nn,chv)
c*****PUTS COMMAND FILE OR INTERACTIVE INPUT INTO TEMPORARY ARRAYS 
c*****FOR SAVING THE FILE IF REQUESTED
      parameter(numt=50)
      character*6 pvar
      character*56 tvar
      integer*2 ix
      character*6 tindex(numt)
      character*80 tinp(numt)
      logical array
      common/lout/tindex,ix,tinp,array
      character*1 ap
      logical chv
      data ap/''''/

        ix=ix+1
        tindex(ix) = pvar
        tinp(ix) = ' '
        if(chv) then
          tinp(ix)=ap//tvar(1:nn)//ap
        else
          tinp(ix)=tvar(1:nn)
        endif

60    return
      end

c***********************************************************************
      subroutine helpm
      character*1 ans,sp
      character*6 var
      character*79 line,blank
      character*87 file1
      data blank/' '/,sp/#20/
c      open(10,file='grid.hlp',status='old')
c
      open(10,file='\pfhelp.tmp',form='formatted',status='old')
      read(10,102) line
  102 format(a)
      close(10)
      n=index(line,' ')
      file1=line(1:n-1)//'grid.hlp'
      open(10,file=file1,form='formatted',status='old')
c
      write(*,'(2a\)')' want general help or a specific variable',
     & '[ g or s] '
      read(*,'(a)') ans

      if(ans.eq.'g'.or.ans.eq.'G') then
        n = 0
   10   read(10,'(a)',end=150) line
        write(*,'(1x,a)') line
        n = n + 1
        if(n.lt.23) go to 10
        write(*,'(a\)') ' more? [y]'
        read(*,'(a)') ans
        if(ans.eq.'n'.or.ans.eq.'N') go to 200
        write(*,'(1x,a)') line
        n = 0
        go to 10
      endif

   20   write(*,'(/a\)') ' enter variable name or ENTER to quit: '
        read(*,'(a)') var
        if(var(1:1).eq.sp) go to 200
        call upcase(var)
        write(*,'(1x,a)') blank
   30   read(10,'(a)',end=100) line
        if(line(2:7).ne.var) go to 30
        write(*,'(1x,a)') line
   40   read(10,'(a)',end=200) line
        if(line.ne.blank) then
          write(*,'(1x,a)') line
          go to 40
        else
          rewind(10)
          go to 20
        endif

  100 write(*,'(a,a6,a)') ' variable: ',var,' not found'
      write(*,'(a\)') ' more help ? [y] '
      read(*,'(a)') ans
      if(ans.ne.'n'.and.ans.ne.'N') then
        rewind(10)
        go to 20
      endif 
      go to 200
  150 pause 'END OF HELP FILE - press ENTER to return to menu'
  200 close(10)
      return
      end
c****************************************************************
      subroutine upcase(var)
C
      character*(*) var
C
      ilen=len(var)
      do 10 i=1,ilen
        k=ichar(var(i:i))
        if(k.ge.97.and.k.le.122) var(i:i)=char(k-32)
   10 continue
      return
      end
