c       GDMERGER
c	Merges two overlapping data sets gridded using the same origin
c	and grid parameters.  The data sets are weighted and summed in
c	the area of overlap.
c
c       written by Jeff Phillips and Tien Grauch
c
c      basic code from program JMRG by Phillips.
c      modified for oddf by Tien Grauch March 1994
c      Grauch added code to construct the large-size matching grids before 
c       merging and plug only in the overlap region Oct 1994
c      modified to restrict overlap areas by Tien Grauch Nov 1994
c
c      parameter nwork is the max rows X cols allowed for overlap region
c
        parameter ( nwork=50000, maxcol=1000 )
        dimension iwork(nwork), work(nwork)
c        equivalence(iwork,work)
        dimension d1(maxcol),d2(maxcol)
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      common /rcspecs/ics1,ics2,ice1,ice2,jrs1,jrs2,jre1,jre2,
     1 icso,iceo,jrso,jreo
	character*50 ifil1,ifil2,ofile
      character state*6, ans*1, id*80,ido*80
        data dval/1.0e+38/
        data dvtest/1.0e+36/
c state is status of temp files upon closing.  change to 'keep'
c  for debugging purposes
c     state='keep'
      state='delete'
cu      call pfinit( 'gdjmrg')
      print 800
800   format(/' Program gdmerger, version 1.0'//
     1' GDMERGER merges two overlapping data sets by weighting'/
     2' then summing in the area of overlap.'/
     3' To accommodate the temporary files generated',
     4' during this program,'/' you will need 3 times the size of',
     5' the output grid + 2 times the'/' size of the overlap between',
     6' grids.'//' Technical contacts - Jeff Phillips and Tien Grauch'/)
1      write(6,100)
100	format(' first grid: '$)
	read(5,105) ifil1
105	format(a)
cu      call gopen(10, ifil1, 'old', 'read', ierr)
        open(10,file=ifil1,form='unformatted',status='old')
cu      if(ierr.ne.0) go to 1
cu      call gh1f4('r', 10, id, nc,nr,xo,dx,yo,dy,ierr)
	read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      print 812, id,xo,yo,nc,nr,dx,dy
812   format(/'title: ',a80/5x,'xo=',g10.4,15x,'yo=',g10.4,/
     1 5x,'nc=',i4,21x,'nr=',i4/5x,'dx=',g10.4,15x,'dy=',g10.4/)
2      write(6,110)
110	format(' second grid: '$)
	read(5,105) ifil2
cu      call gopen(11, ifil2, 'old', 'read', ierr)
        open(11,file=ifil2,form='unformatted',status='old')
cu      if(ierr.ne.0) go to 2
cu      call gh1f4('r', 11, id, nc2,nr2,xo2,dx2,yo2,dy2,ierr)
        read(11) id,pgm,nc2,nr2,nz,xo2,dx2,yo2,dy2
      print 812, id,xo2,yo2,nc2,nr2,dx2,dy2
c
4       print 801
801   format('Merge options:'/5x,'1 - treat both grids with equal',
     1' priority'/5x'2 - restrict overlap region, priority to grid 1',
     2/5x,'3 - restrict overlap region, priority to grid 2'/
     3/10x,'enter option: ',$)
      read*,iopt
      if(iopt.lt.1.or.iopt.gt.3) go to 4
      nband=0
      if(iopt.gt.1) then
        print 802
802     format('Restrict overlap region to how many grid pts? ',$)
        read*,nband
        if(nband.le.0) nband=1
      endif
c
	write(6,130)
130	format(' output grid: '$)
	read(5,105) ofile
      print*,'Enter output title'
      read 106, ido
106   format(a80)
c
c check grid nodes  -- this section inserted by Tien from gdbiharm2.f
      if(dx.ne.dx2.or.dy.ne.dy2) stop 'grid intervals unequal'
      print*
cu      type*,'xo,xo2,dx=',xo,xo2,dx
cu      type*,'yo,yo2,dy=',yo,yo2,dy
      print*,'xo,xo2,dx=',xo,xo2,dx
      print*,'yo,yo2,dy=',yo,yo2,dy
      if(amod((xo-xo2),dx).gt.0.00001) then
        print 814
814   format('grid nodes do not match (roundoff?).'/
     1' Continue anyway? ',$)
        read 105, ans
        if(ans.ne.'y'.and.ans.ne.'Y') stop
      else
         if(amod((yo-yo2),dy).gt.0.00001) then
          print 814
          read 105, ans
          if(ans.ne.'y'.and.ans.ne.'Y') stop
         endif
      endif
c
c find specs of output and specs of overlap work area 
c

      call find_specs(xst,yst,nclap,nrlap,xout,yout,ncout,nrout,
     1 iEW,iNS)
      if(ncout.gt.maxcol) then
         print 815, maxcol
815      format('too many output cols, max is ',i5)
         stop
      endif
      if(nclap*nrlap.gt.nwork) then
        print 816, nwork
816     format('merge-boundary work area too big, max nc*nr =',i8)
        stop
      endif
      if(nband.ge.nclap. or.nband.ge.nrlap) then
        print*,'restriction on overlap region too big'
        stop
      endif
c
c  output grid1 & grid2 data into big grids , units 14,15, resp.
c  non-defined grid nodes are set to dvals
c
cu      call gopen(14, 'jmrg1.tmp', 'new', ' ', ierr)
cu      call gopen(15, 'jmrg2.tmp', 'new', ' ', ierr)
cu      call gh1f4('w', 14, id,ncout,nrout,xout,dx,yout,dy,ierr)
cu      call gh1f4('w', 15, id,ncout,nrout,xout,dx,yout,dy,ierr)
      open(14,form='unformatted',status='scratch')
      write(14)id,pgm,ncout,nrout,nz,xout,dx,yout,dy
      open(15,form='unformatted',status='scratch')
      write(15)id,pgm,ncout,nrout,nz,xout,dx,yout,dy
c construct large grid1
      jr1=0
      do 220 j=1,nrout

      do 211 k=1,maxcol
      d2(k)=dval
211   continue

      jrow=j
      if(jrow.ge.jrs1.and.jrow.le.jre1) then
         jr1=jr1+1
cu         call growf4('r', 10, jr1, d1, nc, ierr)
         read(10) jr1,(d1(i),i=1,nc)
         ic1=0
         do 215 i=1,ncout 
         icol=i
         if(icol.ge.ics1.and.icol.le.ice1) then
            ic1=ic1+1
            d2(icol)=d1(ic1)
         endif
215      continue
      endif

cu      call growf4('w', 14, jrow, d2, ncout, ierr)
      write(14) jrow,(d2(i),i=1,ncout)

220   continue
c
c construct large grid2
      jr2=0
      do 320 j=1,nrout

      do 311 k=1,maxcol
      d1(k)=dval
311   continue

      jrow=j
      if(jrow.ge.jrs2.and.jrow.le.jre2) then
         jr2=jr2+1
cu         call growf4('r', 11, jr2, d2, nc2, ierr)
         read(11) jr2,(d2(i),i=1,nc2)
         ic2=0
         do 315 i=1,ncout 
         icol=i
         if(icol.ge.ics2.and.icol.le.ice2) then
            ic2=ic2+1
            d1(icol)=d2(ic2)
         endif
315      continue
      endif

cu      call growf4('w', 15, jrow, d1, ncout, ierr)
      write(15) jrow,(d1(i),i=1,ncout)

320   continue
cu      call gclose(10, 'keep')
cu      call gclose(11, 'keep')
cu      call gclose(14, 'keep')
cu      call gclose(15, 'keep')
cu      call gopen(14, 'jmrg1.tmp', 'old', ' ', ierr)
cu      call gopen(15, 'jmrg2.tmp', 'old', ' ', ierr)
cu      call gh1f4('r', 14, id,ncout,nrout,xout,dx,yout,dy,ierr)
cu      call gh1f4('r', 15, id,ncout,nrout,xout,dx,yout,dy,ierr)

      close(10)
      close(11)
      rewind(14)
      rewind(15)
      read(14) id
      read(15) id

c  construct weighting grid for overlap area
c
cu      call gopen(12, 'jmrg.tmp', 'new', ' ', ierr)
cu      call gh1f4('w', 12, id,nclap,nrlap,xst,dx,yst,dy,ierr)
      open(12,form='unformatted',status='scratch')
      write(12) id,pgm,nclap,nrlap,nz,xst,dx,yst,dy
        print*
        print *,'...assigning weights of 0, 1, and unknown...'
c
c      weight of 1.0  - data available in grid1 only
c                0.0  - data available in grid2 only
c                dval - data available in both grids
c                0.5  - data available in neither grid     
c
c
c ngrid1 and ngrid2 are temporary flags to tell if grid1, grid2
c  have valid data in them or not.  ialldvals is a flag that
c  checks to see if there is any valid data
c    mncol,mnrow and mxcol,mxrow are values that will be used
c    to restrict the search when readjusting weights (if applicable)
c    they describe approx. where the area of dvals (overlap) are.
c
      mncol=999999
      mnrow=999999
      mxcol=-999999
      mxrow=-999999
      ialldvals=1
      jlap=0
      do 400 j=jrso, jreo
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      read(14) jrow,(d1(i),i=1,ncout)
      read(15) jrow,(d2(i),i=1,ncout)
      jlap=jlap+1

      ilap=0
      nrowpt= (jlap-1) * nclap

      do 350 i=icso, iceo
      ilap=ilap+1
      npt= nrowpt + ilap
      ngrid1=0 
        if(d1(i).lt.dvtest) ngrid1=1 
c
      ngrid2=0
        if(d2(i).lt.dvtest) ngrid2=1
c
c note:
c    ngrid1   ngrid2   wt 
c       0       0      .5
c       1       0       1
c       1       1      dval
c       0       1       0  
c
c assemble the work array with all the weights 
      if(ngrid1.eq.0) then
        if(ngrid2.eq.0) then
           work(npt)=0.5e0
        else
           work(npt)=0.0e0
           ialldvals=0
        endif
      else
        if(ngrid2.eq.0) then
           work(npt)=1.0e0
           ialldvals=0
        else
           work(npt)=dval
           ialldvals=0
           mncol=min(mncol,ilap)
           mxcol=max(mxcol,ilap)
           mnrow=min(mnrow,jlap)
           mxrow=max(mxrow,jlap)
        endif
      endif 
c
c
350   continue
cu      call growf4('w', 12, jlap, work(nrowpt+1), nclap, ierr)
      write(12) jlap,(work(nrowpt+i),i=1,nclap)
400   continue
c
      if(ialldvals.ne.0) stop 'no data in overlap area'
      mncol=mncol-1
      mxcol=mxcol+1
      mnrow=mnrow-1
      mxrow=mxrow+1
      if(mncol.lt.1.or.mncol.gt.nclap) mncol=1
      if(mnrow.lt.1.or.mnrow.gt.nrlap) mnrow=1
      if(mxcol.gt.nclap.or.mxcol.lt.1) mxcol=nclap
      if(mxrow.gt.nrlap.or.mxrow.lt.1) mxrow=nrlap
c
c if needed, readjust weights to reflect overlap restriction
c
      if(iopt.eq.1) go to 440
      print*
      print*,'readjusting weights to restrict overlap region...'
      print*
      if(iopt.eq.2) then
        ghipr=1.0e0
        glopr=0.0e0
      else
        ghipr=0.0e0
        glopr=1.0e0
      endif

c            
c find a value indicating data only in the low priority grid 
c  and replace all dvals within an nband radius with a flag 
c  value of 2.0
c
      ireadjust=0
      do 420 jlap=mnrow,mxrow
      nrowpt=(jlap-1)*nclap

      do 415 ilap=mncol,mxcol
      npt= nrowpt + ilap

      if(work(npt).eq.glopr) then
             do k=1, nband
             if(ilap+k.le.mxcol) then 
                nptk=npt + k
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(ilap-k.ge.mncol) then 
                nptk=npt - k
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(jlap+k.le.mxrow) then 
                nptk= npt + k*nclap
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             if(jlap-k.ge.mnrow) then 
                nptk= npt - k*nclap
                if(work(nptk).eq.dval) work(nptk)=2.0e0
             endif
             end do
      endif
415   continue 
420   continue
c
c  replace all values of 2 with dvals and all dvals with the value
c   of the high priority grid and write out new weights
c
      do 430 jlap=1,nrlap
      nrowpt= (jlap-1)*nclap 
      do 425 ilap=1,nclap
      npt= nrowpt + ilap
      if(work(npt).eq.dval) then
         work(npt)=ghipr
         ireadjust=1
         go to 425
      endif
      if(work(npt).eq.2.0) then
         work(npt)=dval
         ireadjust=1
      endif
425   continue
cu      call growf4('w', 12, jlap, work(nrowpt+1), nclap, ierr)
      write(12) jlap,(work(nrowpt+i),i=1,nclap)
430   continue
      if(ireadjust.eq.0) print*,' no readjustment made'
c
c reopen unit 12 in order to read
c new unit 13 will be jmrg.wgt = plugged weighting grid, ready to use in
c       overlap area.
c
cu440   call gclose(12, 'keep')
cu      call gopen(12, 'jmrg.tmp', 'old', ' ', ierr)
cu      call gh1f4('r', 12, id, nclap, nrlap, xst, dx, yst, dy,ierr)
cu      call gopen(13, 'jmrg.wgt', 'new', 'w', ierr)
cu      call gh1f4('w', 13, id,nclap,nrlap,xst,dx,yst,dy,ierr)
440   rewind(12)
      read(12) id
      open(13,file='jmrg.wgt',form='unformatted',status='unknown')
      write(13) id,pgm,nclap,nrlap,nz,xst,dx,yst,dy
c
        nwrk=nclap*nrlap
          call mcplug(20,12,13,ier,work,iwork,nwrk,nclap,nrlap)
c
	if(ier.eq.1) go to 91
c
cu      call gclose(12,state)
cu      call gclose(13,'keep')
cu      call gopen(13, 'jmrg.wgt', 'old', ' ', ierr)
cu      call gh1f4('r', 13, id,nclap,nrlap,xst,dx,yst,dy,ierr)
      close(12)
      rewind(13)
      read(13) id
c unit 18 is final output grid
c note units 14 & 15 are still open
cu      call gopen(18, ofile, 'new', ' ', ierr)
cu      call gh1f4('w', 18,ido,ncout,nrout,xout,dx,yout,dy,ierr)
      open(18,file='ofile',form='unformatted',status='unknown')
      write(18) ido,pgm,ncout,nrout,nz,xout,dx,yout,dy
      print 812, ido,xout,yout,ncout,nrout,dx,dy
c
c start making the final output grid.
c reuse the array work for output
c

      if(jrso.eq.1) go to 455
      do 450 j= 1,jrso-1
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      read(14) jrow,(d1(i),i=1,ncout)
      read(15) jrow,(d2(i),i=1,ncout)

      if(iNS.eq.1.or.iNS.eq.2) then
         do i= 1, ncout 
         work(i) = d1(i)
         enddo
      else
         do i=1, ncout
         work(i) = d2(i)
         enddo
      endif

cu      call growf4('w', 18, jrow, work, ncout, ierr)
      write(18) jrow,(work(i),i=1,ncout)
450   continue
      
c apply weighting in the overlap region

455   jlap=0
      do 600 j=jrso,jreo
      jrow=j
      jlap=jlap+1
cu      call growf4('r', 13, jlap, work(icso), nclap, ierr)
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      read(13) jlap,(work(icso-1+i),i=1,nclap)
      read(14) jrow,(d1(i),i=1,ncout)
      read(15) jrow,(d2(i),i=1,ncout)

      if(icso.eq.1) go to 480
        if(iEW.eq.1.or.iEW.eq.2) then
           do i=1,icso-1
           work(i)= d1(i)
           end do
        else
           do i=1,icso-1
           work(i)= d2(i)
           end do
        endif
c
c overlap region
c
480     ilap=0
        do 490 i=icso, iceo
        if(d1(i).lt.dvtest) then
         if(d2(i).lt.dvtest) then
            wt0= work(i)
            wt1= 1.e0 - wt0
            work(i) = d1(i)*wt0 + d2(i)*wt1
         else
            work(i)=d1(i)
         endif
        else
         if(d2(i).lt.dvtest) then
            work(i) = d2(i)
         else
            work(i)=dval
         endif
        endif
490     continue
c
c rest of columns if applicable
c
        if(iceo.eq.ncout) go to 560
        if(iEW.eq.0.or.iEW.eq.2) then
           do i=iceo+1, ncout
           work(i)= d1(i)
           end do
        else
           do i=iceo+1, ncout
           work(i)= d2(i)
           end do
        endif

cu560   call growf4('w', 18, jrow, work, ncout, ierr)
560   write(18) jrow,(work(i),i=1,ncout)

600   continue
c
c rest of rows if not in overlap region
c
      if(jreo.eq.nrout) go to 660
      do 650 j= jreo+1, nrout
      jrow=j
cu      call growf4('r', 14, jrow, d1, ncout, ierr)
cu      call growf4('r', 15, jrow, d2, ncout, ierr)
      read(14) jrow,(d1(i),i=1,ncout)
      read(15) jrow,(d2(i),i=1,ncout)

      if(iNS.eq.0.or.iNS.eq.2) then
         do i= 1, ncout 
         work(i) = d1(i)
         end do
      else
         do i= 1, ncout 
         work(i) = d2(i)
         end do
      endif
cu      call growf4('w', 18, jrow, work, ncout, ierr)
      write(18) jrow,(work(i),i=1,ncout)
650   continue

cu660   call gclose(14, 'delete')
cu      call gclose(15, 'delete')
cu      call gclose(13, state)
cu      call gclose(18, 'keep')
660   close(14)
      close(15)
      close(13)
      close(18)
	stop
91	print *,' error in mcplug '
	stop
	end
c***********************************************************************
cc
      subroutine mcplug(nim,idv,jdv,ier,zg,iqd,nwrk,nc,nr)
      dimension n(4)
      dimension zg(nwrk),iqd(nwrk)
c      common /plug1/ zg(100000)
c      common /plug2/ iqd(100000)
c  assuming dimension of wz is maxcol Tien
      dimension wz(1000)
cv      data dv/'37777677777'o/,p2/'m-c ','plug'/
      data dv/1.0e+38/, dvtest/1.0e+36/
c      nwrk=100000
      ier=0
c      read(idv)id,p,nc,nr,nz,xo,dx,yo,dy
      nn=nc*nr
      if(nn.gt.nwrk) go to 99
      do 1 i=1,nn
1        iqd(i)=-1
      ndx=1
      do 2 i=1,nr
      jrow=i
cu      call growf4('r', idv, jrow, zg(ndx), nc, ie)
         call rowio(nc,zg(ndx),-1,idv,jdv,ie)
         if(ie.eq.1) go to 99
2        ndx=ndx+nc
      do 3 i=1,nn
         if(zg(i).ge.1.e29) iqd(i)=0
3        continue
c
c  gridr and curvmn are from the minc program
c
20    call gridr(nc,nr,zg,wz,ier)
c      go to 5
      na=0
      do 21 i=1,nn-nc
      if(iqd(i).eq.0) then
        if(zg(i).gt.0.49.and.zg(i).lt.0.51) then
c          iqd(i)=-1
          na=na+1
        else if (zg(i).ne.zg(i+1)) then
          temp=0.5*(zg(i)+zg(i+1))
          if(temp.gt.0.49.and.temp.lt.0.51) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else if (zg(i).ne.zg(i+nc)) then
          temp=0.5*(zg(i)+zg(i+nc))
          if(temp.gt.0.49.and.temp.lt.0.51) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else
          zg(i)=dv
        endif
      endif
21    continue
      print *,na,' unknown weights replaced by 0.5'
c      go to 5
      call gridr(nc,nr,zg,wz,ier)
      na=0
      do 22 i=1,nn-nc
      if(iqd(i).eq.0) then
        if((zg(i).gt.0.24.and.zg(i).lt.0.26).or.(zg(i).gt.0.74.and.
     1  zg(i).lt.0.76)) then
c          iqd(i)=-1
          na=na+1
        else if (zg(i).ne.zg(i+1)) then
          temp=0.5*(zg(i)+zg(i+1))
          if((temp.gt.0.24.and.temp.lt.0.26).or.(temp.gt.0.74.and.
     1    temp.lt.0.76)) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else if (zg(i).ne.zg(i+nc)) then
          temp=0.5*(zg(i)+zg(i+nc))
          if((temp.gt.0.24.and.temp.lt.0.26).or.(temp.gt.0.74.and.
     1    temp.lt.0.76)) then
            zg(i)=temp
c            iqd(i)=-1
            na=na+1
          else
            zg(i)=dv
          endif
        else
          zg(i)=dv
        endif
      endif
22    continue
      print *,na,' unknown weights replaced by 0.25 or 0.75'
      call gridrold(nc,nr,zg,wz,ier)
      if(nim.eq.0) go to 5
23      eps=0.
      call curvmn(zg,iqd,wz,nc,nr,eps,nim,a,b,ni)
c5     write(jdv)id,p2,nc,nr,nz,xo,dx,yo,dy
5      ndx=1
      do 4 i=1,nr
        jrow=i
         call rowio(nc,zg(ndx),0,idv,jdv,ie)
cu         call growf4('w', jdv, jrow, zg(ndx), nc, ie)
4        ndx=ndx+nc
      return
99    ier=1
      return
      end
c*********************************************************************
      subroutine curvmn(zg,iqd,b,nc,nr,epsmx,nim,eps1,dn1,
     1ni)
c
c  Applies minimum curvature equations to the first nc*nr elements of 
c  array zg.  Array iqd contains nc*nr elements which indicate for each 
c  mesh location the quadrant where a data value is located.  An iqd 
c  value of zero indicates no data and -1 locks the present mesh value.
c  Array b should contain 6*nc*nr elements used for weighting when iqd 
c  is 1 to 4, in the case where iqd is only 0 or -1, b can be of length 
c  one.  The over-relaxation parameter w increases as the system 
c  converges until 1.7 is reached.
c
      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
c first row
c
      if(iqd(1))2,1,1
1     zg(1)=(( (2.*(zg(2)+zg(nc+1))-zg(nc+nc+1)-zg(3))*.5 )-
     & 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)-
     & 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)-
     &      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)-
     &   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 )-
     &   zg(nc))*w+zg(nc)
c
c second row
c
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)-
     &   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))-
     &   2.*zg(j1+1)-zg(jm+1)-zg(j1-1)-zg(i+2)-zg(j1+nc))*
     &   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))-
     &      2.*(zg(j1-1)+zg(j1+1))-zg(jm-1)-zg(jm+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)-
     &   zg(jm-1)-zg(j1+1)-zg(i-2)-zg(j1+nc))*5.5555556e-2 )-
     &   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)-
     &   zg(i-2)-zg(j1-1)-zg(j1+nc))*.16666667 )-zg(i))*w+zg(i)
c
c rows 3 to nr-2
c
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)
     &     -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)
     &      -2.*(zg(j1+1)+zg(jm+1))-zg(j1-1)-zg(j1+nc)-zg(i+2)-
     &      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)+zg(jm-1)+zg(j1-1))-zg(j1+nc)-zg(jm-nc)-zg(i-2)-
     &          zg(i+2))*.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(j1+1)+b2*zg(i+1)+b3*zg(jm)+b4*zg(jm-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(jm-1)+b2*zg(i-1)+b3*zg(j1)+b4*zg(j1+1)
32          t=.25*(zg(j1+nc)+zg(i-2)+zg(jm-nc)+zg(i+2))
     &        +.5*(zg(j1-1)+zg(jm-1)+zg(jm+1)+zg(j1+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)+
     &      zg(jm-1))-zg(jm+1)-zg(jm-nc)-zg(i-2)-
     &      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)-zg(j1-1)-zg(j1+nc))*.14285714 )-zg(i))*w+zg(i)
39       continue
c
c row nr-1
c
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)-
     & 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))-
     & 2.*zg(jm+1)-zg(jm-1)-zg(jm-nc)-zg(i+2)-
     & 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)-
     &    2.*(zg(jm-1)+zg(jm+1))-zg(j1-1)-zg(i-2)-
     &    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)-
     & zg(j1-1)-zg(i-2)-zg(jm-nc)-zg(jm+1))*5.5555556e-2 )-
     & 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)-
     & zg(i-2)-zg(j1-1))*.16666667 )-zg(i))*w+zg(i)
c
c last row
c
50    i=i+1
      if(iqd(i))52,51,51
51    jm=i-nc
      zg(i)=(( (2.*(zg(i+1)+zg(jm))-zg(i+1)-zg(jm-nc))*.5 )-
     & 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)-
     & 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)-
     &    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)-
     & 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 )-
     & 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 gridr(nc,nr,zg,wz,ier)
      dimension zg(nc,nr),w(1)
10    call mediplug(nc,nr,zg,nd)
      if(nd.gt.0) go to 10
      return
      end
c     mediplug.for - dos version
c     J.Phillips, USGS, Denver, Colorado
c     changes from vms version are indicated by cv
c
c     fills holes using the median of surrounding values
c************************************************************************
c
      subroutine mediplug(nc,nr,zg,nd)
      dimension zg(1)
      dimension a(1000), b(1000, 3), c(1000), id(14)
      character ifile*50, id*56
      data ddval / 1.0e+36 /
c    1 write(*, 200)
c  200 format(14h enter ifile: ,$)
c      read(5, 100) ifile
c  100 format(a50)
cv    open(10,file=ifile,status='old',form='unformatted')
c      open(10, file=ifile, status='old', form='unformatted', err=1)
c      write(*, 201)
c  201 format(14h enter ofile: ,$)
c      read(5, 100) ifile
cv    open(11,file=ifile,status='new',form='unformatted')
c      open(11, file=ifile, status='unknown', form='unformatted')
c      read(10) id, pgm, nc, nr, nz, xo, dx, yo, dy
c      write(11) id, pgm, nc, nr, nz, xo, dx, yo, dy
      do 100 i=1,nc
      b(i,3)=zg(i)
      b(i,2)=zg(i+nc)
  100 continue
c      read(10) dlt, (b(i,3),i = 1, nc)
c      write(11) dlt, (b(i,3),i = 1, nc)
c      read(10) dlt, (b(i,2),i = 1, nc)
      nd = 0
      do 14 ii = 3, nr
      do 120 i=1,nc
      b(i,1)=zg(i+(ii-1)*nc)
  120 continue
c        read(10) dlt, (b(i,1),i = 1, nc)
        c(1) = b(1,2)
        c(nc) = b(nc,2)
        do 12 i = 2, nc - 1
          c(i) = b(i,2)
          if (c(i) .gt. ddval) then
            j = 0
            do 11 k = 1, 3
              do 10 l = i - 1, i + 1
                if (b(l,k) .lt. ddval) then
                  j = j + 1
                  a(j) = b(l,k)
                end if
10            continue
11          continue
            if (j .gt. 0) call median(j, a)
            if (j .lt. 1) then
              nd = nd + 1
            else
              c(i) = a(1)
            end if
          end if
12      continue
        do 130 i=1,nc
        zg(i+(ii-2)*nc) = c(i)
130     continue
c        write(11) dlt, (c(i),i = 1, nc)
        do 13 i = 1, nc
          b(i,3) = b(i,2)
          b(i,2) = b(i,1)
13      continue
14    continue
c      write(11) dlt, (b(i,1),i = 1, nc)
c      write(*, 101) nd
  101 format(i10,13h dvals remain)
      return
      end
c
c***************************************************************************
      subroutine median(n, a)
      dimension a(n)
   10 if (n .le. 2) return 
      ax = a(1)
      an = a(1)
      iax = 1
      ian = 1
      do 11 i = 2, n
        if (a(i) .gt. ax) then
          ax = a(i)
          iax = i
        end if
        if (a(i) .lt. an) then
          an = a(i)
          ian = i
        end if
11    continue
      if (an .eq. ax) then
        a(1) = an
        return 
      end if
      j = 0
      do 20 i = 1, n
      temp = a(i)
      if ((i .eq. iax) .or. (i .eq. ian)) goto 20
      j = j + 1
      a(j) = temp
   20 continue
      n = j
      goto 10
      end
c*****************************************************************************
      subroutine gridrold(nc,nr,zg,wz,ier)
      dimension zg(1),wz(1)
      data dval/1.0e+38/,ddval/1.0e+36/
      nn=nc*nr
      i=1
      if(zg(1).ge.ddval) then
        do 10 i=2,nn
        if(zg(i).lt.ddval) go to 20
10      continue
        stop 'empty grid'
20      do 30 j=1,i-1
        zg(j)=zg(i)
30      continue
      endif
      do 70 j=i,nn
      if(zg(j).lt.ddval) go to 70
      if(j.eq.nn) go to 51
      do 50 k=j+1,nn
      if(zg(k).gt.ddval) go to 50
      slope=(zg(k)-zg(j-1))/float(k-j+1)
      do 40 ii=j,k-1
      zg(ii)=zg(j-1)+slope*float(ii-j+1)
40    continue
      go to 70
50    continue
51    continue
      do 60 ii=j,nn
      zg(ii)=zg(j-1)
60    continue
      return
70    continue
      end
c
c *************************************************************************
      subroutine find_specs(xst,yst,nclap,nrlap,xout,yout,ncout,nrout,
     1 iEW,iNS)
c subroutine to find the overlap area of two grids from the grid specs
c in order to minimize the size of grid to work on.
c   adapted from gdbiharm2.f by Tien
c
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      common /rcspecs/ics1,ics2,ice1,ice2,jrs1,jrs2,jre1,jre2,
     1 icso,iceo,jrso,jreo
c
c ----OVERLAP SPECS--------
c
c find general area where grids touch/overlap
      xx= xo + float(nc - 1)*dx
      yx= yo + float(nr - 1)*dy
      xx2= xo2 + float(nc2 - 1)*dx
      yx2= yo2 + float(nr2 - 1)*dy
c
c CASES
c  directional codes below say where grid1 is in relation to grid2
c    iEW=0  grid1 on E
c    iEW=1  grid1 on W
c    iEW=2  grid1 on E and W
c    iEW=3  grid2 on E and W
c    iNS=0  grid1 on N
c    iNS=1  grid1 on S
c    iNS=2  grid1 on N and S
c    iNS=3  grid2 on N and S
c
c CASE    iEW or iNS   overlap?
c  A          1           yes
c  B          0           yes
c  C          1            no    --stop program
c  D          0            no    --stop program
c  E          2           yes
c  F          3           yes
c
c Determine starting, ending x & y values of grid overlap
c based on the grid specifications
c
c
      dx2=2.0e0*dx
      dy2=2.0e0*dy
c
c EAST-WEST CASES
c take care of round off problems in the difference
      xodiff=(xo-xo2)
      if(abs(xodiff).le.0.00001) xodiff=0.0
      xxdiff=(xx-xx2)
      if(abs(xxdiff).le.0.00001) xxdiff=0.0
      xxo2diff=(xx-xo2)
      if(abs(xxo2diff).le.0.00001) xxo2diff=0.0
      xx2odiff=(xx2-xo)
      if(abs(xx2odiff).le.0.00001) xx2odiff=0.0
c  Case A
      if(xodiff.le.0.0.and.xxdiff.le.0.0.and.xxo2diff.ge.0.0) then
         iEW=1
         xst=amax1(xo,xo2-dx2)
         xend=amin1(xx2,xx+dx2)
         go to 7
      endif
c  Case B
      if(xodiff.ge.0.0.and.xxdiff.ge.0.0.and.xx2odiff.ge.0.0) then
         iEW=0
         xst=amax1(xo2,xo-dx2)
         xend=amin1(xx,xx2+dx2)
         go to 7
      endif
c  Case E
      if(xodiff.lt.0.0.and.xxdiff.gt.0.0.and.xxo2diff.gt.0.0) then
         iEW=2
         xst=amax1(xo,xo2-dx2)
         xend=amin1(xx,xx2+dx2)
          go to 7
      endif
c  Case F
      if(xodiff.gt.0.0.and.xxdiff.lt.0.0.and.xx2odiff.gt.0.0) then
         iEW=3
         xst=amax1(xo2,xo-dx2)
         xend=amin1(xx2,xx+dx2)
         go to 7
      endif
c  Case C and D
      print 810
810   format(' Grids have no overlap in x.  Please extend',
     1 ' grid2 to overlap with grid1')
      stop
c
c NORTH-SOUTH CASES
c
c take care of round off problems in the difference
7     yodiff=(yo-yo2)
      if(abs(yodiff).le.0.00001) yodiff=0.0
      yxdiff=(yx-yx2)
      if(abs(yxdiff).le.0.00001) yxdiff=0.0
      yxo2diff=(yx-yo2)
      if(abs(yxo2diff).le.0.00001) yxo2diff=0.0
      yx2odiff=(yx2-yo)
      if(abs(yx2odiff).le.0.00001) yx2odiff=0.0
c  Case A
      if(yodiff.le.0.0.and.yxdiff.le.0.0.and.yxo2diff.ge.0.0) then
         iNS=1
         yst=amax1(yo,yo2-dy2)
         yend=amin1(yx2,yx+dy2)
         go to 8
      endif
c  Case B
      if(yodiff.ge.0.0.and.yxdiff.ge.0.0.and.yx2odiff.ge.0.0) then
         iNS=0
         yst=amax1(yo2,yo-dy2)
         yend=amin1(yx,yx2+dy2)
        go to 8
      endif
c  Case E
      if(yodiff.lt.0.0.and.yxdiff.gt.0.0.and.yxo2diff.gt.0.0) then
         iNS=2
         yst=amax1(yo,yo2-dy2)
         yend=amin1(yx,yx2+dy2)
         go to 8
      endif
c  Case F
      if(yodiff.gt.0.0.and.yxdiff.lt.0.0.and.yx2odiff.gt.0.0) then
         iNS=3
         yst=amax1(yo2,yo-dy2)
         yend=amin1(yx2,yx+dy2)
         go to 8
      endif
c  Case C and D
      print 811
811   format(' Grids have no overlap in y.  Please extend',
     1 ' grid2 to overlap with grid1')
      stop
c  find no. of cols and rows of overlap region 
8       nclap=nint( (xend - xst) / dx + 1.0 )
        nrlap=nint( (yend - yst) / dy + 1.0 )
c
c
c--------OUTPUT SPECS-------
c
c find specs of output merged (large) grid
c
        xout=amin1(xo,xo2)
        xxout=amax1(xx,xx2)
        yout=amin1(yo,yo2)
        yxout=amax1(yx,yx2)
        ncout= nint((xxout-xout)/dx + 1.e0)
        nrout= nint((yxout-yout)/dy + 1.e0)
c
c  ---FIND WHERE THESE AREAS ARE RELATIVE TO OUTPUT GRID--
c
c find the rows and cols in output grid assoc. with 
c beginning and ending grid1 and grid2
c    ics1,ice1= col start,end of grid 1 in large grid
c    jrs1,jre1= row start,end of grid 1 in large grid
c    ics2,ice2= col start,end of grid 2 in large grid
c    jrs2,jre2= row start,end of grid 2 in large grid
c
      call rcfind( xout, dx, xo, ics1)
      call rcfind( xout, dx, xx, ice1)
      call rcfind( yout, dy, yo, jrs1)
      call rcfind( yout, dy, yx, jre1)
      call rcfind( xout, dx, xo2, ics2)
      call rcfind( xout, dx, xx2, ice2)
      call rcfind( yout, dy, yo2, jrs2)
      call rcfind( yout, dy, yx2, jre2)
c
c find cols and rows in output grid where overlap area begins and ends
        call rcfind(xout, dx, xst, icso)
        call rcfind(xout, dx, xend, iceo)
        call rcfind(yout, dy, yst, jrso)
        call rcfind(yout, dy, yend, jreo)
c
c if any of overlap is outside the output grid, readjust overlap area
c
c also check overlap region to see if it could be at edge of merged grid
c if so, there may be no control for plugging
c
        if(icso.lt.1) then
          if(jrso.le.1.and.jreo.ge.nrout) call warning
          xst=xst + dx * float( 1 - icso)
          icso=1
        endif
        if(iceo.gt.ncout) then
          if(jrso.le.1.and.jreo.ge.nrout) call warning
          nclap=nclap - (iceo - ncout)
          iceo=ncout
        endif
        if(jrso.lt.1) then
          yst=yst + dy * float( 1 - jrso)
          jrso=1
        endif
        if(jreo.gt.nrout) then
           nrlap=nrlap - (jreo - nrout)
           jreo=nrout
        endif
c
        print*
        print*,'xo,yo,nc,nr of merge-boundary work area=',
     1xst,yst,nclap,nrlap
        print*,'xo,yo,nc,nr of final output grid=',
     1xout,yout,ncout,nrout
c       type*,'beg,end cols & rows of grid1:',ics1,ice1,jrs1,jre1
c       type*,'beg,end cols & rows of grid2:',ics2,ice2,jrs2,jre2
c       type*,'beg,end cols & rows of overlap:',icso,iceo,jrso,jreo
        print*
        return
        end
 
c***************************************************************************
      subroutine rcfind(xyo, del, f, nf)
 
      ff= (f - xyo)/del + 1.e0
      nf= nint( ff )
 
      return
      end
c****************************************************************************
      subroutine warning
c warning message subroutine
c
        character*1 ans
        write (6,132)
132     format(/'WARNING overlap area may extend to edge of grid.'/
     1  5x,'Use the subset ("su") option of GDUTILITY to add a',
     2 ' border',/,5x,'of dvals to the input grids if so.'/)
          write(6,133)
133       format(' do you want to continue?: '$)
          read(5,105) ans
105   format(a)
          if(ans.eq.'n'.or.ans.eq.'N') stop 'OK'
        return
        end
c****************************************************************************
      subroutine rowio(n,z[huge],iop,idev,jdev,iend)
c
c  WHERE IOP<0 READ; IOP=0 WRITE; IOP>0 READ&WRITE
c
      dimension z(n)
      iend=0
      if(iop)1,2,1
1     read(idev,end=10) xo,z
      if(iop)9,9,2
2     write(jdev) xo,z
9     return
10    iend=1
      return
      end
