

c             program gmerge
c    Merge two grids without interpolation, no restriction on placement
c  or grid dimensions.  The increments must be the same and grid nodes
c  must register.
c
c    Intersection is handled in two modes:
c  1) First grid has priority and second grid only overwrites dvals.
c  2) Second grid has priority and overwrites the contributions of the 
c     first grid regardless of the contents.
c
c  Mike Webring, USGS
c    version 1.1, update 7/90.
c    changed logic in xfer to match overwriting scheme documented above.
c    separated common areas, removed enddo's, and tab continuations.
c
	common /grdspc/ nc,  nr,  xo,  dx,  yo,  dy,
     1                  nc2, nr2, xo2, dx2, yo2, dy2,
     1                  nc3, nr3, xo3, dx3, yo3, dy3,
     1                  isw, jsw, ksw, dval, iproj, cmdeg, bldeg
	common /grdtit/ id,    id2,    id3,    pgm,   pgm2,   pgm3
	character       id*56, id2*56, id3*56, pgm*8, pgm2*8, pgm3*8

	parameter       ( nwork=15000 )
	dimension       work(nwork)
	character*80    f1, f2, f3

1	format(a)
	dval = 1.7e37
	nhdr = 1

	isw  = 10
	jsw  = 11
	ksw  = 12

	print *, ' version 1.1'
10	print *, ' first input grid :'
	read(5,1) f1
	open( unit=isw, file=f1, status='old', form='unformatted',
     1     mode='read', iostat =ios )
	if ( ios .ne. 0 ) then
	  print *, ' error opening', f1
	  go to 10
	endif

20	print *, ' second input grid :'
	read(5,1) f2
	open ( unit=jsw, file=f2, status='old', form='unformatted',
     1     mode='read', iostat=ios )
	if ( ios .ne. 0 ) then
	  print *, ' error opening ', f2
	  go to 20
	endif

	call setf3

	print *, ' output columns and rows =', nc3, nr3
	if ( nc + nc2 + nc3  .gt.  nwork ) then
	  print *, ' ncol + ncol2 + ncol_out >', nwork
	  stop
	endif
	if ( nc3 .gt. 200 + nc + nc2 ) then
	  print *, ' the gap between 1 and 2 is greater than 200 columns'
	  print *, ' do you want to continue ? [y]'
	  ians = noyes()
	  if ( ians .eq. 0  .or.  ians .eq. -2 ) stop
	endif
	if ( nr3 .gt. 200 + nr + nr2 ) then
	  print *, ' the gap between 1 and 2 is greater than 200 rows'
	  print *, ' do you want to continue ? [y]'
	  ians = noyes()
	  if ( ians .eq. 0  .or.  ians .eq. -2 ) stop
	endif

	id3  = id
	pgm3 = pgm

	print *, ' merged output grid :'
	read(5,1) f3
	open ( unit=ksw, file=f3, status='new', form='unformatted' )

	write( ksw ) id3, pgm3, nc3,nr3,nhdr, xo3,dx3, yo3,dy3, 
     1            iproj, cmdeg, bldeg

	n2 =  1 + nc
	n3 = n2 + nc2
	call rowctl( work, work(n2), work(n3) )

	stop
	end
c*****************************************************************************
	subroutine setf3

	common /grdspc/ nc,  nr,  xo,  dx,  yo,  dy,
     1                  nc2, nr2, xo2, dx2, yo2, dy2,
     1                  nc3, nr3, xo3, dx3, yo3, dy3,
     1                  isw, jsw, ksw, dval, iproj, cmdeg, bldeg
	common /grdtit/ id,    id2,    id3,    pgm,   pgm2,   pgm3
	character       id*56, id2*56, id3*56, pgm*8, pgm2*8, pgm3*8

	call grdhdr( isw, id, pgm,   nc,nr,nz,   xo,dx,   yo,dy,
     1               iproj, cmdeg, bldeg )
	call grdhdr( jsw, id2, pgm2, nc2,nr2,nz, xo2,dx2, yo2,dy2,
     1               iproj2, cm2, bl2 )

	if ( dx .eq. 0.0 ) stop ' dx = 0'
	if ( dy .eq. 0.0 ) stop ' dy = 0'

c  check increments
	testdx = abs( ( dx - dx2 ) / dx )
	testdy = abs( ( dy - dy2 ) / dy )
	if ( testdx .gt. 0.001 ) stop ' dx.ne.dx2'
	if ( testdy .gt. 0.001 ) stop ' dy.ne.dy2'

c  flip coordinate system to right cartesian
	flipx = 1.0
	flipy = 1.0
	if ( dx .lt. 0.0 ) then
	  xo    = -xo
	  xo2   = -xo2
	  dx    = -dx
	  flipx = -1.0
	endif
	if ( dy .lt. 0.0 ) then
	  yo    = -yo
	  yo2   = -yo2
	  dy    = -dy
	  flipy = -1.0
	endif

c  check x registration
	istop = 0
	gu    = abs ( ( xo2 - xo ) / dx )
	xt    = gu - aint( gu )
	if ( xt .le. 0.9  .and.  xt .ge. 0.1 ) then
	  print *, ' xo, xo2, grid unit offset', xo, xo2, xt
	  print *, ' x misregistered more than .1 unit'
	  istop = 1
	endif

c  check y registration
	gu = abs ( ( yo2 - yo ) / dy )
	yt = gu - aint( gu )
	if ( yt .le. 0.9  .and.  yt .ge. 0.1 ) then
	  print *, ' yo, yo2, grid unit offset', yo, yo2, yt
	  print *, ' y misregistered more than .1 unit'
	  istop = 1
	endif
	if ( istop .eq. 1 ) stop

c  find min/max x
	xm1  = xo  + dx  * float( nc  - 1 )
	xm2  = xo2 + dx2 * float( nc2 - 1 )
	xmin = xo
	if ( xo2 .lt. xo ) xmin = xo2
	xmax = xm1
	if ( xm2 .gt. xm1 ) xmax = xm2

c  find min/max y
	ym1  = yo  +  dy * float( nr  - 1 )
	ym2  = yo2 + dy2 * float( nr2 - 1 )
	ymin = yo
	if ( yo2 .lt. yo ) ymin = yo2
	ymax = ym1
	if ( ym2 .gt. ym1 ) ymax = ym2

c  output grid specs
	nc3 = 1 + int( ( ( xmax - xmin ) / dx ) + 0.1 )
	nr3 = 1 + int( ( ( ymax - ymin ) / dy ) + 0.1 )
	xo3 = flipx * xmin
	yo3 = flipy * ymin
	dx3 = flipx * dx
	dy3 = flipy * dy

	return
	end
c*****************************************************************************
	subroutine rowctl( r, r2, r3 )

	common /grdspc/ nc,  nr,  xo,  dx,  yo,  dy,
     1                  nc2, nr2, xo2, dx2, yo2, dy2,
     1                  nc3, nr3, xo3, dx3, yo3, dy3,
     1                  isw, jsw, ksw, dval, iproj, cmdeg, bldeg
	common /grdtit/ id,    id2,    id3,    pgm,   pgm2,   pgm3
	character       id*56, id2*56, id3*56, pgm*8, pgm2*8, pgm3*8

	dimension  r(nc), r2(nc2), r3(nc3)
	logical    read1, read2, firstg, ovrwrt

	fltbig = 1.0e35

	print *,' does the first grid have priority ? [y]'
	ovrwrt = .false.
	if ( noyes() .eq. 0 ) ovrwrt = .true.
	print *, ' second grid overwriting =', ovrwrt

	is1    = 1 + int( ( ( xo  - xo3 ) / dx ) + 0.1 )
	is2    = 1 + int( ( ( xo2 - xo3 ) / dx ) + 0.1 )

	dyerr  = abs( 0.1 * dy )
	read1  = .true.
	read2  = .true.
	irow1  = 0
	irow2  = 0
	ieof1  = 0
	ieof2  = 0

	do 100 iout = 1, nr3

	yout   = yo3  + dy3 * float( iout - 1 )
	ylower = yout - dyerr
	yupper = yout + dyerr

	do 20 i = 1, nc3
	  r3(i) = dval
20	continue

c  input row 1
	if ( read1   .and.   irow1 .lt. nr ) then
	  irow1 = irow1 + 1
	  call rowio( nc, r, -1, isw, isw, ieof1 )
	  if ( ieof1 .eq. 0 ) then
	    y1    = yo + dy * float( irow1 - 1 )
	    read1 = .false.
	    do 30 i = 1, nc
	      if ( r(i) .gt. fltbig ) r(i) = dval
30	    continue
	  else
	    irow1 = irow1 - 1
	    y1    = dval
	    print *, ' EOF file1 after row', irow1
	  endif
	endif

c  input row 2
	if ( read2   .and.   irow2 .lt. nr2 ) then
	  irow2 = irow2 + 1
	  call rowio( nc2, r2, -1, jsw, jsw, ieof2 )
	  if ( ieof2 .eq. 0 ) then
	    y2    = yo2 + dy * float( irow2 - 1 )
	    read2 = .false.
	    do 40 i = 1, nc2
	     if ( r2(i) .gt. fltbig ) r2(i) = dval
40	    continue
	  else
	    irow2 = irow2 - 1
	    y2    = dval
	    print *, ' EOF file2 after row', irow2
	  endif
	endif

c  insert row 1
	if( y1 .gt. ylower   .and.   y1 .lt. yupper  ) then
	  firstg = .true.
          call xfer( is1, nc, r, nc3, r3, firstg, ovrwrt )
	  read1 = .true.
	endif

c  insert row 2
	if( y2 .gt. ylower   .and.   y2 .lt. yupper ) then
	  firstg = .false.
          call xfer( is2, nc2, r2, nc3, r3, firstg, ovrwrt )
	  read2 = .true.
	endif

c  output row
	write( ksw ) yout, r3

100	continue

	if ( irow1 .ne. nr ) print *,
     1  'error: file 1, input vs nrow =', irow1, nr
	if ( irow2 .ne. nr2 ) print *,
     1 'error: file 2, input vs nrow =', irow1, nr

	return
	end
c******************************************************************************
        subroutine xfer( istart, nc, r, nc3, r3,
     1                   firstg, ovrwrt )

	dimension  r(nc), r3(nc3)
	logical    firstg, ovrwrt

	fltbig = 1.0e30

	iend   = istart - 1 + nc
	if ( istart .le. 0   .or.
     1       iend   .gt. nc3  ) stop ' xfer: index error'
	j     = istart

	if ( firstg ) then

	  do 20 i = 1, nc
	    r3(j) = r(i)
	    j     = j + 1
20	  continue
	
	else

	  do 30 i = 1, nc
	    if ( ovrwrt ) then
	      r3(j) = r(i)
	    else
	      if ( r3(j) .gt. fltbig ) r3(j) = r(i)
	    endif
	    j = j + 1
30	  continue

	endif

	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
c*****************************************************************************
	subroutine grdhdr( iunit, id, pgm, nc,nr,nz, xo,dx,yo,dy,
     1                        iproj, cmerid, baslat )
c  return info contained in the grid header record.
c  projection info is optional.
	character id*56, pgm*8
	read( iunit, err=9 )  id, pgm, nc,nr,nz, xo,dx,yo,dy,
     1                     iproj, cmerid, baslat
	return
9	rewind iunit
	read( iunit, err=99 ) id, pgm, nc,nr,nz, xo,dx,yo,dy
	iproj  = 0
	cmerid = 999.
	baslat = 999.
	return
99	print *, ' cannot read grid header'
	stop
	end
c****************************************************************************
	function noyes()
c  returns -2 for ctl-Z,   -1 for no change,
c           0 for no,      1 for yes

	character inp*80, blank*80

	blank = ' '
	icnt  = 0
	noyes = -1

1	inp   = ' '
	read ( 5, 2, end=9, err=9 ) inp
2	format( a80 )
	if ( inp .eq. blank ) return

	do 3 n = 1, 4
	  if ( inp(n:n) .ne. ' ' ) then
	    if ( inp(n:n).eq.'y' .or. inp(n:n).eq.'Y' ) noyes = 1
	    if ( inp(n:n).eq.'n' .or. inp(n:n).eq.'N' ) noyes = 0
	    if ( noyes .gt. -1 ) return
	    if ( ichar( inp(n:n) ) .le. 31 ) go to 9
	    go to 5
	  endif
3	continue

5	if ( icnt .lt. 3 ) then
	  icnt = icnt + 1
	  write(6,6)
6	  format(' no,yes, or <cr>:'$)
	  go to 1
	endif
	write(6,*) ' count exceeded...taking default'
	noyes = -1
	return

c  abort
9	noyes = -2
	return
	end
