c*********************************************************************
cu      subroutine curvmn(zg,iqd,b,nc,nr,epsmx,nim,eps1,dn1,
      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
c************************************************************************
cu      subroutine gridr(nc,nr,zg,wz,ier)
      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************************************************************************
c
cu      subroutine mediplug(nc,nr,zg,nd)
      subroutine mediplug(nc,nr,zg[huge],nd)
      dimension zg(1)
      dimension a(1000), b(1000, 3), c(1000)
c      dimension id(14)
c      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*****************************************************************************
cu      subroutine gridrold(nc,nr,zg,wz,ier)
      subroutine gridrold(nc,nr,zg[huge])
c      dimension zg(1),wz(1)
      dimension zg(1)
c      data dval/1.0e+38/
      data 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****************************************************************************
c      subroutine rowio(n,z[huge],iop,idev,jdev,iend)
cc
cc  WHERE IOP<0 READ; IOP=0 WRITE; IOP>0 READ&WRITE
cc
c      dimension z(n)
c      iend=0
c      if(iop)1,2,1
c1     read(idev,end=10) xo,z
c      if(iop)9,9,2
c2     write(jdev) xo,z
c9     return
c10    iend=1
c      return
c      end
