c                program jigsaw

c  trim multiple grids to polygonal boundaries defined in the command file.
c  the resulting grids are the minimum size required to contain the remaining
c  data, and no-data values (variable 'dv') are inserted where neccesary to
c  conform to the polygon.  once the grids are trimmed to a set of consistent
c  boundaries, they can be merged using a simple program.
c
c namelist parameters
c
c xvert	array of x coordinates in data units (kilometers, degrees or
c	whatever units in the grid header record)
c
c yvert	array of y coordinates 
c
c nvert	number of elements in x/yvert, limit is 200
c
c vfile	optional file containing x,y coordinate pairs.
c	either vfile or xvert,yvert,nvert can be specified.
c	vfile is read in free-field format with one coordinate pair per line.
c
c cutbac  distance in data units (km, feet, etc) matching those of the grid
c	header record.  a positive cutbac inserts dvals inside the specified
c	boundary while a negative does the opposite to allow overlap of grids.
c
c ngrid	number of grid specs following the namelist
c
c
c example 
c
c 3  1--------2---------------3       Three separate grids referenced to the
c    |         \         2    |       same coordinate system.  
c 2  |    1     7-----8-------9       
c    |           \       3    |       
c 1  4------------5-----------6
c
c    1  2  3  4  5  6  7  8  9
c
c &parms
c xvert= 1, 4, 9.5, 1, 5.5, 9.5, 4.5, 6.5, 9.5
c yvert= 3, 3,   3, 1,   1,    1,  2,   2,   2
c nvert=9,cutbac=1.,ngrid=3
c &
c test1.grd   << file name of the input grid
c 5           << number of vertices in the next line
c 1 2 7 5 4   << list of vertices defining a polygon in clockwise order
c test2.grd
c 5
c 2 3 9 8 7
c test3.grd
c 5
c 5 7 8 9 6
c
c  author: mike webring, usgs
c
	common /jigsaw/ z(2000)
	dimension list(200),xvert(200),yvert(200),xv(200),yv(200)
	character*56 ifile,ofile,vfile,id
	character ext*3,p*8
	logical inside
	common /parms/ vfile,xvert,yvert,cutbac,ngrid,nvert
	data dv/0.1701412e+39/, ext/'jig'/
c  dv is a large floating point number approximately equal to 10**37.
c
	vfile=' '
	write (*,1)  
1	format(' enter command filename :'$)
	read(5,37) ifile
	open(unit=9,file=ifile,status='old',form='formatted',mode='read')
        call namemc(9)
	if(vfile(1:1).ne.' ') then
	open(unit=12,file=vfile,form='formatted',status='old',mode='read')
	nvert=0
5	nvert=nvert+1
	if(nvert.gt.200) stop ' vertex list > 200 entries'
	read(12,*,end=6) xvert(nvert),yvert(nvert)
	go to 5
6	nvert=nvert-1
	close(12)
	endif
c
	do 200 k=1,ngrid
	read(9,37) ifile
37	format(a56)
	read(9,*) ncrnr
	read(9,*) (list(i),i=1,ncrnr)       
	if(ncrnr.le.2) go to 200
	open(unit=13,file=ifile,status='old',form='unformatted',
     1 mode='read')
	n=index(ifile,'.')
	ofile=ifile(1:n)//ext
	open(unit=14,file=ofile,status='unknown',form='unformatted')
	do 50 i=1,ncrnr
	xv(i)=xvert(list(i))
50	yv(i)=yvert(list(i))
c
	twoc=2.0*cutbac
	if (cutbac.eq.0.0) go to 29
	a1=areap(ncrnr,xv,yv)
	do 21 i=1,ncrnr
	xv(i)=xv(i)+cutbac
	a2=areap(ncrnr,xv,yv)
	if(cutbac*(a1-a2).ge.0.0) go to 21
	xv(i)=xv(i)-twoc
	a2=areap(ncrnr,xv,yv)
21	a1=a2
	do 22 i=1,ncrnr
	yv(i)=yv(i)+cutbac
	a2=areap(ncrnr,xv,yv)
	if(cutbac*(a1-a2).ge.0.0) go to 22
	yv(i)=yv(i)-twoc
	a2=areap(ncrnr,xv,yv)
22	a1=a2
c
29	read(13) id,p,nc,nr,nz,xo,dx,yo,dy
	xmx=xv(1)
	xmn=xv(1)
	ymx=yv(1)
	ymn=yv(1)
	do 100 i=2,ncrnr
	if(xv(i).gt.xmx) xmx=xv(i)
	if(xv(i).lt.xmn) xmn=xv(i)
	if(yv(i).gt.ymx) ymx=yv(i)
	if(yv(i).lt.ymn) ymn=yv(i)
100	continue
	xmn=amax1(xo,xmn)
	x2=dx*float(nc-1)+xo
	xmx=amin1(x2,xmx)
	ymn=amax1(yo,ymn)
	y2=dy*float(nr-1)+yo
	ymx=amin1(y2,ymx)
c
	ix=int((xmn-xo)/dx)
	xo2=xo+dx*float(ix)
	iy=int((ymn-yo)/dy)
	yo2=yo+dy*float(iy)
	nc2=int((xmx-xmn)/dx)+1
	nr2=int((ymx-ymn)/dy)+1
	if(ix.lt.0) stop 111
	if(iy.lt.0) stop 222
	if(nc2.gt.nc) stop 333
	if(nr2.gt.nr) stop 444
	if(nc2.lt.0) nc2=1
	if(nr2.le.0) nr2=1
	write(14) id,p,nc2,nr2,nz,xo2,dx,yo2,dy
	ixs=ix+1
	ixe=ix+nc2
	iys=iy+1
	iye=iy+nr2
c
	if(iys.gt.1) then
	do 110 j=1,iys-1
110	call rowio(nc,z,-1,13,14,ie)
	endif
c
	yp=yo2
	do 40 j=iys,iye
	call rowio(nc,z,-1,13,14,ie)
	xp=xo2
	do 30 i=ixs,ixe
	if(z(i).gt.1.e29) z(i)=dv
	if(.not.inside(ncrnr,xv,yv,xp,yp)) z(i)=dv
30	xp=xp+dx
	call rowio(nc2,z(ixs),0,13,14,ie)
40	yp=yp+dy
	close(13)
	close(14)
200	continue
c
	close(9)
	stop
	end
c******************************************************************************
	function areap(n,x,y)
	dimension x(n),y(n)
	areap=0.0
	if(n.lt.3) return
	xm=x(1)
	do 1 i=2,n
	if(x(i).gt.xm) xm=x(i)
1	continue
	x1=xm-x(1)
	do 2 i=1,n
	i1=i+1
	if(i1.gt.n) i1=1
	x2=xm-x(i1)
	areap=areap+0.5*(y(i1)-y(i))*(x1+x2)
2	x1=x2
	return
	end
c******************************************************************************
	logical function inside(n,x,z,xp,zp)
c   is xp,zp inside polygon defined by x,z ?
	dimension x(1),z(1)
	logical z1,z2
	inside=.false.
	do 5 i=1,n
	i1=i+1
	if(i1.gt.n) i1=1
	dz=z(i1)-z(i)
	if(dz.eq.0.0) go to 5
	z1 = zp.le.z(i)
	z2 = zp.gt.z(i1)
	if((z1.and.z2) .or. (.not.z1 .and. .not.z2)) go to 3
	go to 5
c   points on boundary are outside with condition .lt.
c   for inclusion change to .le.
3	rslope=(x(i1)-x(i))/dz
	d=(xp-x(i))-(zp-z(i))*rslope
	if(d.lt.0.0) inside = .not.inside
5	continue
	return
	end
c******************************************************************************
	subroutine rowio(n,z,iop,idev,jdev,iend)
c  read iop<0, write iop=0, r&w iop=1
	dimension z(n)
	iend=0
	if(iop)1,2,1
1	read(idev,end=10) y,z
	if(iop)9,9,2
2	write(jdev) y,z
9	return
10	iend=1
	return
	end

      subroutine namemc(icmd)
c
c     namelist simulator for pc's with no namelist statement
c     this subroutine is program independent but it
c     calls subroutine check which is program dependent
c     The program does not handle the repeat indicator *
c     nor does it handle subscripted array variables or
c     substrings of character variables.
c     It also does not handle complex or logical variables.
c     nn = number of characters in value name(tvar)
c     chv = character variable indicator(logical)
c     inum = an array index number
c     nvar = a number returned by subroutine check to
c     determine whether a variable is an array
c     numa = a number returned by subroutine check
c     which is compared against nvar
c     r.godson,usgs,denver,co., 11/87
c
      parameter(ivar=80)
      character var(ivar),comma,apos,blank,equal,amp,dollar
      character*6 pvar
      character*56 tvar
      logical chv,vset
      data comma/','/,apos/#27/,blank/' '/,equal/'='/
      data amp/'&'/,dollar/'$'/,pvar/' '/,tvar/' '/
c
c     get namelist start name(e.g.&parms)
c
    5 read(icmd,1000,end=910) var
      do 7 i=1,ivar
      if(var(i).ne.amp.and.var(i).ne.dollar) go to 7
      do 6 m=i+1,ivar
      if(var(m).eq.blank) go to 15
    6 continue
    7 continue
      go to 5
c
c     start processing variables
c
   10 read(icmd,1000,end=900) var
 1000 format(80a1)
      m=1
   15 do 20 i=m,ivar
      if(var(i).eq.amp.or.var(i).eq.dollar) go to 900
      if(var(i).ne.blank.and.var(i).ne.comma) go to 30
   20 continue
      go to 10
   30 mm=0
c
c     check for continuation of array values on succeding lines
c
      item=ichar(var(i))
      if(item.eq.43.or.item.eq.45.or.item
     & .eq.46.or.(item.gt.47.and.item.lt.58).or.var(i).eq.apos) then
      k=i
      go to 95
      endif
c
c     get program variable name
c
   35 inum=1
      pvar=blank
      vset=.false.
      do 40 j=i,ivar
      if(var(j).eq.blank.or.var(j).eq.equal) go to 50
      mm=mm+1
      pvar(mm:mm)=var(j)
   40 continue
c
c     get variable value
c
     

   50 do 80 k=j+1,ivar
      if(var(k).eq.blank.or.var(k).eq.equal) go to 80
      if(var(k).eq.apos) then
c
c     character variable
      chv=.true.
      vset=.true.
      nn=0
      do 60 l=k+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
   60 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      else
c
c     not a character variable
      chv=.false.
      nn=0
      do 70 l=k,ivar
      if(var(l).eq.blank.or.var(l).eq.comma
     & .or.var(l).eq.amp.or.var(l).eq.dollar) go to 90
c
c     check for non numeric character
c
      item=ichar(var(l))
      if((item.lt.48.and.(item.ne.46
     & .and.item.ne.45.and.item.ne.43)) .or. (item .gt. 57 .and.
     & (item .ne. 69 .and. item .ne. 101 .and. item.ne.68
     &  .and. item.ne.100))) then
      write(*,*) 'missing first apostropy or a non numeric',
     & ' character in namelist variable ',pvar
      stop
      endif
      nn=nn+1
      tvar(nn:nn)=var(l)
   70 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no delimiter'
      stop
      endif
   80 continue
      write(*,*) 'error in namelist variable ',pvar,' no  value'
      stop
c
c     call program dependent subroutine to assign values
c
   90 call check(pvar,tvar,nn,chv,nvar,numa,inum)
      tvar=blank
      k=l+1
c
c     check to see if variable is an array
c     the variable numa used is program dependent
c
      if(nvar.lt.numa) go to 110
   95 if(.not.chv) then
c
c     array variable
c
      nn=0
      do 100 l=k,ivar
      if(var(l).eq.blank.and.nn.eq.0) go to 100
      if(var(l-1).eq.blank.and.var(l).eq.comma.and.nn.eq.0) go to 100
      if((var(l).eq.amp.or.var(l).eq.dollar).and.nn.eq.0) go to 900
c
c     check for consecutive commas
c
      if(var(l).eq.comma.and.vset) then
      vset=.false.
      inum=inum+1
      go to 90
      else if(var(l).eq.comma) then
      inum=inum + 1
      go to 100
      endif
c
      if(var(l).eq.blank
     & .or.var(l).eq.amp.or.var(l).eq.dollar) then
      inum=inum+1
      go to 90
      endif
      if(ichar(var(l)).gt.57.and.nn.eq.0) go to 120
      nn=nn+1
      tvar(nn:nn)=var(l)
      vset=.true.
  100 continue
      else
c
c     character array variable
c
      do 105 j=k,ivar
      if(var(j).eq.blank) go to 105
      if(var(j).eq.amp.or.var(j).eq.dollar) go to 900
      if(ichar(var(j)).gt.57) then
      l=j
      go to 120
      endif
c
c     check for consecutive commas
c
      if(var(j).eq.comma.and.vset) then
      vset=.false.
      go to 105
      else if(var(j).eq.comma) then
      inum=inum + 1
      go to 105
      endif
      if(var(j).eq.apos) then
      nn=0
      inum=inum+1
      vset=.true.
      do 103 l=j+1,ivar
      if(var(l).eq.apos) go to 90
      nn=nn+1
      tvar(nn:nn)=var(l)
  103 continue
      write(*,*) 'error in namelist variable ',pvar,
     & ' no ending apostrophy'
      stop
      endif
  105 continue
c
c     end of array values
c
      endif
  110 m=l+1
      if(var(l).eq.amp.or.var(l).eq.dollar) go to 900
      go to 130
  120 m=l
  130 if(m.lt.81) then
      go to 15
      else
      go to 10
      endif
  900 return
  910 write(*,*) ' error in namelist;no beginning & or $'
      stop
      end
        subroutine check(pvar,tvar,nn,chv,nvar,numa,inum)
c
c       assigns values to proper variables
c       variables are passed to program, jigsaw, through common blocks
c       numr=position in array var where real variables start
c       numa=position in array var where arrays start
c       nnvar=number of variables in program jigsaw
c
        parameter (nnvar=6,numr=3)
        character*6 pvar,var(nnvar)
        character*56 tvar,kvar,cfmt
        logical chv
        character*56 vfile
        dimension xvert(200),yvert(200)
        common/parms/vfile,xvert,yvert,cutbac,ngrid,nvert
        data var/'ngrid','nvert','cutbac','xvert','yvert','vfile'/
        numa=4
        do 190 i=1,nnvar
        if (pvar .ne. var(i)) goto 190
        if (.not. chv) then
c
c       noncharacter value
c       right justify the number in variable kvar
c
        m=57-nn
        im=m-1
        kvar(m:56)=tvar(1:nn)
        if (i .lt. numr) then
c
c       integer value
c
        write(cfmt,50) im,nn
50      format('(',i2,'x,i',i2,')')
        read(kvar,cfmt) jvar
        else
c
c       real value
c
        write(cfmt,60) im,nn
60      format('(',i2,'x,g',i2,'.0)')
        read(kvar,cfmt) xvar
        endif
        endif
        goto (101,102,103,104,105,106),i
101     ngrid=jvar
        goto 200
102     nvert=jvar
        goto 200
103     cutbac=xvar
        goto 200
104     xvert(inum)=xvar
        goto 200
105     yvert(inum)=xvar
        goto 200
106     vfile=tvar(1:nn)
        goto 200
190     continue
        write(*,*)' error in namelist - ',pvar,' variable not included'
        stop
200     nvar=i
        return
        end
