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
c
        parameter ( nwork=50000 )
        dimension iwork(nwork), work(nwork)
c        equivalence(iwork,work)
        dimension id(14),pgm(2),d1(1000),d2(1000),d1sav(1000)
	character*50 ifil1,ifil2,ofile
        character*1 ans
cv        data dval/'37777677777'o/
        data dval/1.0e+38/,ddval/1.0e+30/
	write(6,100)
100	format(' first grid: '$)
	read(5,105) ifil1
105     format(a)
        open(10,file=ifil1,form='unformatted',status='old')
	write(6,110)
110	format(' second grid: '$)
	read(5,105) ifil2
        open(11,file=ifil2,form='unformatted',status='old')
	write(6,130)
130	format(' output grid: '$)
        read(5,105) ofile
        open(14,file=ofile,form='unformatted',status='new')
c        open(12,file=ofile,form='unformatted',status='new')
c        close(12,status='delete')
c        open(12,form='unformatted',status='scratch')
c        open(12,file='jmrg.tmp',form='unformatted',status='unknown')
c
	read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
	read(11) id,pgm,nc1,nr1,nz,xo1,dx1,yo1,dy1
c        if(nc.ne.nc1.or.nr.ne.nr1.or.xo.ne.xo1.or.yo.ne.yo1.or.dx.ne.dx1
c     1       .or.dy.ne.dy1) goto 90
        if(nc.ne.nc1.or.nr.ne.nr1) go to 90
c        write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
        write(14) id,pgm,nc,nr,nz,xo,dx,yo,dy
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
c        if(nband.le.0) nband=1
        if(nband.eq.0) nband=1
        if(iopt.eq.2) then
          ghipr=1.0e0
          glopr=0.0e0
        else
          ghipr=0.0e0
          glopr=1.0e0
        endif
      endif
c
        write(6,131)
131     format(' do you want to save the weighting grid?: '$)
        read(5,105) ans
        if(ans.eq.'y'.or.ans.eq.'Y') then
          ans='y'
          open(15,file='jmrg.wgt',form='unformatted',status='unknown')
          write(15) id,pgm,nc,nr,nz,xo,dx,yo,dy
          print *,'weighting grid saved in file "jmrg.wgt"'
        endif
        ncnr=nc*(nr+1)
        nbnd=ncnr/nwork+1
        print*,'Processing in ',nbnd,' bands'
        jbnd=nr/nbnd+1
        jb=1
        do 80 kk=1,nbnd
        je=min(nr,jb+jbnd-1)
        print*,'Band ',kk,': Rows ',jb,' to ',je
        open(12,form='unformatted',status='scratch')
        nr1=je-jb+1
        if(jb.eq.1) then
          write(12) id,pgm,nc,nr1,nz,xo,dx,yo,dy
        else
          write(12) id,pgm,nc,nr1+1,nz,xo,dx,yo,dy
          write(12) dlt,(d1sav(i),i=1,nc)
        endif
c
c  construct weighting grid
c
        print *,'assigning weights of 0, 1, and unknown'
        iflag=0
c        do 20 j=1,nr
        rewind(10)
        rewind(11)
        read(10) id
        read(11) id
        if(jb.ne.1) then
          do 19 j=1,jb-1
          read(10) dlt
          read(11) dlt
19        continue
        endif
        do 20 j=jb,je
	read(10) dlt,(d1(i),i=1,nc)
	read(11) dlt,(d2(i),i=1,nc)
	do 10 i=1,nc
        if (d1(i).lt.ddval) then
                if (d2(i).lt.ddval) then
c			overlap area
			d1(i) = dval
                        if(i.eq.1.or.i.eq.nc.or.j.eq.1.or.j.eq.nr) then
                         iflag=1
                        endif
			go to 10
		else
c			first grid data area
			d1(i) = 1.0
			go to 10
		end if
	end if
        if (d2(i).lt.ddval) then
c		second grid data area
		d1(i) = 0.0
		go to 10
	end if
c	no data area
	d1(i) = 0.5
10	continue
	write(12) dlt,(d1(i),i=1,nc)
20	continue
        if(iflag.eq.1) then
          print *,' '
          print *,'WARNING overlap area extends to edge of grid.'
          print *,'        Use the extract ("x") option of UTILITY to'
          print *,'        add a border of dvals to the input grids.'
          print *,' '
          write(6,133)
133       format(' do you want to continue?: '$)
          read(5,105) ans
          if(ans.eq.'n'.or.ans.eq.'N') stop
        endif
c        stop
	rewind 10
	rewind 11
	rewind 12
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*
      read(12) id,pgm,nc,nr2,nz,xo,dx,yo,dy
c      jbb=max(1,jb-nband)
c      if(jb.ne.1.and.jbb.gt.1) then
c        do j=1,jbb-1
c        read(12) dlt
c        enddo
c      endif
c      jeb=min(nr,je+nband)
      k=1
c      do j=jbb,jeb
      do j=1,nr2
c      if(k.gt.nwork) stop 'nband is too large'
      read(12) dlt,(work(i),i=k,k+nc-1)
      k=k+nc
      enddo
      close(12)
      open(12,form='unformatted',status='scratch')
      write(12) id,pgm,nc,nr2,nz,xo,dx,yo,dy
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
c      do 420 j=jbb,jeb
      do 420 j=1,nr2
      nrowpt=(j-1)*nc

      do 415 i=1,nc
      npt= nrowpt + i

      if(work(npt).eq.glopr) then
      if(nband.gt.0) then
             do k=1, nband
             if(i+k.le.nc) then
                nptk=npt + k
                if(work(nptk).ge.ddval) work(nptk)=2.0e0
             endif
             if(i-k.ge.1) then
                nptk=npt - k
                if(work(nptk).ge.ddval) work(nptk)=2.0e0
             endif
             if(j+k.le.nr2) then
                nptk= npt + k*nc
                if(work(nptk).ge.ddval) work(nptk)=2.0e0
             endif
             if(j-k.ge.1) then
                nptk= npt - k*nc
                if(work(nptk).ge.ddval) work(nptk)=2.0e0
             endif
             enddo
      else
             do k=1, -nband
             if(i+k.le.nc) then
               nptk=npt + k
               if(work(nptk).ge.ddval) work(npt)=4.0e0
             endif
             if(i-k.ge.1) then
               nptk=npt - k
               if(work(nptk).ge.ddval) work(npt)=4.0e0
             endif
             if(j+k.le.nr2) then
               nptk= npt + k*nc
               if(work(nptk).ge.ddval) work(npt)=4.0e0
             endif
             if(j-k.ge.1) then
                nptk= npt - k*nc
                if(work(nptk).ge.ddval) work(npt)=4.0e0
             endif
             end do
      endif
      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
      k=1
c      do 430 j=jbb,jeb
      do 430 j=1,nr2
      nrowpt= (j-1)*nc
      do 425 i=1,nc
      npt= nrowpt + i
      if(work(npt).ge.ddval) then
         work(npt)=ghipr
         ireadjust=1
         go to 425
      endif
      if(work(npt).eq.2.0e0) then
         work(npt)=dval
         ireadjust=1
      endif
425   continue
      write(12) dlt,(work(i),i=k,k+nc-1)
      k=k+nc
430   continue
c      if(ireadjust.eq.0) print*,' no readjustment made'
      rewind(12)
440   continue
c
          open(13,form='unformatted',status='scratch')
c
c  plug weighting grid
c
        call mcplug(20,12,13,ier,work,iwork,nwork)
	if(ier.eq.1) go to 91
c        stop
	rewind 13
c
c  multiply first grid by weights
c
	close(12)
        open(12,form='unformatted',status='scratch')
        write(12) id,pgm,nc,nr1,nz,xo,dx,yo,dy
	read(10) id
	read(13) id
        if(jb.ne.1) then
          do 39 j=1,jb-1
          read(10) dlt
39        continue
        endif
c        do 40 j=1,nr
        do 40 j=jb,je
	read(10) dlt,(d1(i),i=1,nc)
	read(13) dlt,(d2(i),i=1,nc)
        if(ans.eq.'y') then
          write(15) dlt,(d2(i),i=1,nc)
        endif
	do 30 i=1,nc
        if(d2(i).gt.1.0e0) then
          d1(i)=dval
c          print*,j,i
        endif
        if(d1(i).lt.ddval) d1(i)=d1(i)*d2(i)
30	continue
40	write(12) dlt,(d1(i),i=1,nc)
        do 41 i=1,nc
        d1sav(i)=d2(i)
41      continue
c
c  weight second grid and add to weighted first grid
c
	rewind 12
	rewind 13
        rewind 10
c        close(10)
c        open(10,file=ofile,form='unformatted',status='new')
c        write(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
        read(10) id
	read(11) id
	read(12) id
	read(13) id
        if(jb.ne.1) then
          do 69 j=1,jb-1
          read(11) dlt
69        continue
        endif
c        do 70 j=1,nr
        do 70 j=jb,je
	read(11) dlt,(d1(i),i=1,nc)
	read(13) dlt,(d2(i),i=1,nc)
c	weight second grid
	do 50 i=1,nc
        if(d2(i).gt.1.0e0) then
          d1(i)=dval
c          print*,j,i
        endif
        if(d1(i).lt.ddval) d1(i)=d1(i)*(1.-d2(i))
50	continue
c	read weighted first grid
	read(12) dlt,(d2(i),i=1,nc)
c	add them
	do 60 i=1,nc
	if(d1(i).lt.dval) then
		if(d2(i).lt.dval) d1(i) = d1(i) + d2(i)
		go to 60
	end if
	if(d2(i).lt.dval) then
		d1(i) = d2(i)
		go to 60
	end if
60	continue
70      write(14) dlt,(d1(i),i=1,nc)
        close(12)
        close(13)
        jb=je+1
80      continue
	stop
90	print *,' grids must have identical dimensions and origins '
	stop
91	print *,' error in mcplug '
	stop
	end
      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
      subroutine mcplug(nim,idv,jdv,ier,zg[huge],iqd[huge],nwrk)
      dimension id(14),p(2),p2(2)
      dimension zg(nwrk),iqd(nwrk)
c      common /plug1/ zg(100000)
c      common /plug2/ iqd(100000)
      dimension wz(2000)
cv      data dv/'37777677777'o/,p2/'m-c ','plug'/
      data dv/1.0e+38/,p2/'m-c ','plug'/
c      nwrk=100000
      ier=0
      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
         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  curvmn is from the minc program
c
20    call gridr(nc,nr,zg)
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)
      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)
      if(nim.eq.0) go to 5
      eps=0.
      call curvmn(zg,iqd,wz,nc,nr,eps,nim,a,b,ni)
5     write(jdv)id,p2,nc,nr,nz,xo,dx,yo,dy
      ndx=1
      do 4 i=1,nr
         do 6 j=1,nc
         if(zg(ndx+j-1).eq.4.0e0) go to 6
         if(zg(ndx+j-1).gt.1.0) zg(ndx+j-1)=1.0
         if(zg(ndx+j-1).lt.0.0) zg(ndx+j-1)=0.0
6        continue
         call rowio(nc,zg(ndx),0,idv,jdv,ie)
4        ndx=ndx+nc
      return
99    ier=1
      return
      end
      subroutine curvmn(zg[huge],iqd[huge],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
      subroutine gridr(nc,nr,zg[huge])
      dimension zg(nc,nr)
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
      subroutine mediplug(nc,nr,zg[huge],nd)
      dimension zg(1)
      dimension a(1000), b(1000, 3), c(1000)
c      character ifile*50
      data ddval / 1.0e+30 /
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
      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
      subroutine gridrold(nc,nr,zg[huge])
      dimension zg(1)
      data ddval/1.0e+30/
      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

