c biharm.for - version 3.0
c              Jeff Phillips, November 30, 1995
c
c  version 3.0 has the following changes from version 2.0
c   - outputs statistics for edges of both grids as well as for the
c     overlap region
c   - computes median and (crude) mode of shifts, as well as the mean
c   - outputs the difference grid, as well as the overlap grids
c   - outputs the edge shifts as xyz files rather than as grids
c  version 2.0 has the following changes from version 1.0:
c   - revamp of code to make program faster and better able to handle dvals
c   - additional output of grids of the overlap region between grids.  As 
c      this program only outputs shift values at the survey boundary area,
c      it may be instructive to compare the overlap region in entirety 
c      using the program gdlscale, which gives least squares constant
c      (shift) and slope of the two overlap grids.
c  Please report any abnormalities to Tien Grauch, tien@musette.cr.usgs.gov
c    July 5, 1994
c
c  Description:
c  This program uses the biharmonic operator to check for the constant
c  shift between two grids (in preparation for merging) per Misac
c  Nabighian's method.
c  Input: grid1 is properly trimmed to define the
c  the survey boundaries desired.  Grid2 does not need to be trimmed.
c  Output is a grid of constants determined from the biharmonic operator
c  for the touching region of the two grids.
c
c The biharmonic operator uses the following 13 grid pts:
c
c                          u(1)
c
c                u(2)      u(3)       u(4)
c      
c       u(5)     u(6)      u(7)       u(8)       u(9)
c
c                u(10)     u(11)      u(12)
c
c                          u(13)
c
c The biharmonic operator is
c  20*u(7) - 8*(u(3)+u(8)+u(11)+u(6)) + 2(u(2)+u(4)+u(12)+u(10) + u(1)+
c       u(9)+u(13)+u(5) = bih
c To check the datum shift required to match grid2 to grid1, put
c u(7) on grid1 next to the boundary and the others split between the two
c grids to be merged.  If there is no datum shift between the two grids,
c  bih=0.
c  If there is a shift, shft= bih/(sum of weights for pts lying in grid 2)
c
c  also outputs grids of common data area (or bombs out if no overlap):
c   (1) grid1 for overlap area   grid1lap.grd
c   (2) grid2 for overlap area   grid2lap.grd
c
c   (3) map of data coverage     maplap.gd (temporary)
c  
      parameter (ncmax=3000)
      dimension g1(ncmax), g2(ncmax)
      dimension s1mode(100), s2mode(100), dmode(100)
cu      character igrid1*50, igrid2*50, ofile*50, title1*80
      character igrid1*50, igrid2*50, ofile*50, title1*56
cu      character title2*80, title3*80, id*8, pfile*80, a1*1
      character title2*56, title3*56, pfile*80, a1*1
      common /data/ g1lap(ncmax,5),g2lap(ncmax,5),gmap(ncmax,5)
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      common /bih/ title3,dval,ncout,xoout,yoout,minc,maxc,minr,maxr
      common /sums/ nsum1,nsum2
      dval=1.0e+38
      dvtest=1.0e+30

      call pfinit( 'biharm' )
      print 890
890   format(/' biharm v. 3.0',//
     1 '    This program uses a biharmonic operator to estimate the'
     2 /' shifts along the boundary between 2 grids, as per'
     3 /' Nabighian''s Laplacian method.  It also estimates average'
     4 /' shifts in the area of overlap.  At least one input grid must'
     5 /' be trimmed to the desired merge boundary.')
      print 891
891   format(
     1 '    Output of the program includes xyz files of shift values at'
     2/' the survey boundaries, and three grids that contain the grid1'
     3/' and grid2 values for the overlap region and their difference.'
     4/)
cu      print 892
cu892   format(' Please report any abnormalities to Tien Grauch, ',
cu     1' tien@musette.cr.usgs.gov'//)

1     print*,'Enter name of the first grid'
      read 880, igrid1
880   format(a50)
      call gopen( 10, igrid1, 'old', 'read', ierror )
      if (ierror. ne. 0) then
        print 881
881   format('Error opening grid.  Try again')
        go to 1
      endif
cu      call gh1f4( 'r', 10, title1, nc, nr, xo, dx, yo, dy, ierr)
      call gheader('r',10,title1,nc,nr,xo,dx,yo,dy,ierr)
cu      print 882,xo,yo,dx,nc,nr
882   format(' xo = ',g14.7,2x,'yo = ',g14.7,'   dx,dy = ',g14.7,
     1 /' nc = ',i9,5x,'nr = ',i9/)
5     print*,'Enter name of the second grid'
      read 880, igrid2
      call gopen(11, igrid2, 'old', 'read', ierror)
      if (ierror. ne. 0 ) then
        print 881
        go to 5
      endif
cu      call gh1f4('r',11,title2, nc2, nr2, xo2, dx2, yo2, dy2, ierr2)
      call gheader('r',11,title2,nc2,nr2,xo2,dx2,yo2,dy2,ierr2)
cu      print 882, xo2,yo2,dx2,nc2,nr2

c  check that grid1 and grid2 grid nodes match

      if(ierr. ne. 0. or. ierr2. ne. 0) stop 'problem reading grids'

      if(dx.ne.dx2.or.dy.ne.dy2) stop 'grid intervals unequal' 
      if(dx.ne.dy) stop 'dx must = dy'
cu      print*,'checking grids'
cu      type*,'xo,xo2,dx=',xo,xo2,dx
cu      type*,'yo,yo2,dy=',yo,yo2,dy
cu      print*,'xo,xo2,dx=',xo,xo2,dx,amode((xo-xo2),dx)
cu      print*,'yo,yo2,dy=',yo,yo2,dy,amode((yo-yo2),dy)
cu      if(abs(amod((xo-xo2),dx)).gt.0.000 1) then
      if(abs(amode((xo-xo2),dx)).gt.0.0001) then
          print*,'xo,xo2,dx=',xo,xo2,dx,amode((xo-xo2),dx)
          print*, 'grid nodes must match.  Continue anyway?'
          read 898, a1
898   format(a1)
          if(a1.ne.'y'.and.a1.ne.'Y') stop
      else
cu         if(abs(amode((yo-yo2),dy)).gt.0.00001) then
         if(abs(amod((yo-yo2),dy)).gt.0.0001) then
          print*,'yo,yo2,dy=',yo,yo2,dy,amode((yo-yo2),dy)
          print*, 'grid nodes must match.  Continue anyway?'
          read 898, a1
          if(a1.ne.'y'.and.a1.ne.'Y') stop
         endif
      endif

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  (**NOTE: the directional cases are not used, but may need them
c    some day to speed up the operation???)
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, adding or subtracting
c 2 grid points, to account for needing to reach out two
c grid points on either side of u(7) in the computation of 
c the biharmonic operator
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
c     if(xo.le.xo2. and. xx.le.xx2. and. xx.ge.xo2) then
         iEW=1
         xstart=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
c     if(xo.ge.xo2. and. xx.ge.xx2. and. xx2.ge.xo) then
         iEW=0
         xstart=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
c     if(xo.lt.xo2. and. xx.gt.xx2. and. xx.gt.xo2) then
         iEW=2
         xstart=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
c     if(xo.gt.xo2. and. xx.lt.xx2. and. xx2.gt.xo) then
         iEW=3
         xstart=amax1(xo2,xo-dx2)
         xend=amin1(xx2,xx+dx2)
         go to 7
      endif
c Case G

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
c     if(yo.le.yo2. and. yx.le.yx2. and. yx.ge.yo2) then
         iNS=1
         ystart=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
c     if(yo.ge.yo2. and. yx.ge.yx2. and. yx2.ge.yo) then
         iNS=0
         ystart=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
c     if(yo.lt.yo2. and. yx.gt.yx2. and. yx.gt.yo2) then
         iNS=2
         ystart=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
c     if(yo.gt.yo2. and. yx.lt.yx2. and. yx2.gt.yo) then
         iNS=3
         ystart=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 the rows and cols in grid1 assoc. with xstart,xend,ystart,yend 
c    ics1,ice1= col start,end in grid 1
c    jrs1,jre1= row start,end in grid 1
c    ics2,ice2= col start,end in grid 2
c    jrs2,jre2= row start,end in grid 2
c  okay if these cols, rows are outside grids for now
8     call rcfind( xo, dx, xstart, ics1)
      call rcfind( xo, dx, xend, ice1)
      call rcfind( yo, dy, ystart, jrs1)
      call rcfind( yo, dy, yend, jre1)
c find the equivalent cols and rows in grid2. 
      call rcfind( xo2, dx, xstart, ics2)
      call rcfind( xo2, dx, xend, ice2)
      call rcfind( yo2, dy, ystart, jrs2)
      call rcfind( yo2, dy, yend, jre2)
c note that range of cols or rows can't both be outside 
c grid1 and grid2 
c     if(nrange(ics1,1).ne.0) print*,'ics1 outside grid1'
c     if(nrange(ice1,1).ne.0) print*,'ice1 outside grid1'
c     if(nrange(ics2,3).ne.0) print*,'ics2 outside grid2'
c     if(nrange(ice2,3).ne.0) print*,'ice2 outside grid2'
c     if(nrange(jrs1,2).ne.0) print*,'jrs1 outside grid1'
c     if(nrange(jre1,2).ne.0) print*,'jre1 outside grid1'
c     if(nrange(jrs2,4).ne.0) print*,'jrs2 outside grid2'
c     if(nrange(jre2,4).ne.0) print*,'jre2 outside grid2'
c
c  initialize output grids.  
        ncout=ice1 - ics1 + 1
        nrout=jre1 - jrs1 + 1
c double check col & row range
        ncout2=ice2 - ics2 + 1
        nrout2=jre2 - jrs2 + 1
        if(ncout2.ne.ncout.or.nrout2.ne.nrout) then
            stop 'bug in program--no. cols or rows not same'
        endif
        xoout= xstart
        yoout= ystart
cu        print*
cu        print*,'xo,yo,nc,nr of output=',xstart,ystart,ncout,nrout
cu        print*
cu10      print*,'Enter output grid of shifts'
cu        read 880,ofile
cu        print*,'Enter title'
cu        read 803,title3
cu803     format(a80)
cu        call gopen(15, ofile, 'new', 'write', ierr)
        open(17,form='unformatted',status='scratch')
cu        if(ierr.ne.0) then
cu           print*,'Error opening grid.  Try again'
cu           go to 10
cu        endif
cu        call gh1f4('w',15,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
c
      print*,'Shift values will be saved in files "grid1edg.xyz" and "gr
     1id2edg.xyz"'
      print 813
813   format(/' In addition, 3 grids will be output:'/
     1 '   (1) data of grid1 in overlap area -->  grid1lap.grd'/
     2 '   (2) data of grid2 in overlap area -->  grid2lap.grd'/
     3 '   (3) difference (grid1 - grid2)    -->  g1mg2lap.grd'/
     4 /' Computing...'/)
cu        call gopen(12, 'grid1lap.gd', 'new', 'write', ierr)
cu        call gopen(13, 'grid2lap.gd', 'new', 'write', ierr)
cu        call gopen(14, 'maplap.gd', 'new', 'write', ierr)
        call gopen(12,'grid1lap.grd','unknown','readwrite',ierr)
        call gopen(13,'grid2lap.grd','unknown','readwrite',ierr)
        call gopen(14,' ','scratch','readwrite',ierr)
        call gopen(19,' ','scratch','readwrite',ierr)
        call gopen(16,'g1mg2lap.grd','unknown','write',ierr)
      title3='data from grid1 common to grid2'
cu        call gh1f4('w',12,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      call gheader('w',12,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      title3='data from grid2 common to grid1'
cu        call gh1f4('w',13,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      call gheader('w',13,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
cu        call gh1f4('w',14,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      title3='difference (grid1 - grid2) in overlap area'
      call gheader('w',16,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
c
c set up map of data coverage (maplap.grd) and write other overlap grd's
c   0=no data in either grid1 nor grid2
c   1=data available from grid1 
c   2=data available from grid2 but not grid1
c
c as assemble the overlap grids, keep track of min, max rows & cols
c where data occurs 
c
      nbig=1000000
      minc=nbig
      maxc=-nbig
      minr=nbig
      maxr=-nbig
      jstart=0
      ndif=0
      jr1=jrs1-1
      jr2=jrs2-1
c     type*,'Beginning jr1,jr2=',jr1,jr2

c if the rows start outside of either grid1 or grid2, need to put
c dvals into the associated overlap grid and the proper mlap value
c into the data map grid
      if(jrs1.le.0.and.jrs2.le.0) then
         stop 'starting row outside both grids--bug in program'
      endif
      if(jre1.gt.nr.and.jre2.gt.nr2) then
         stop 'ending rows outside both grids--bug in program'
      endif
c
c case where starting row is outside grid1 only
      if(jrs1.le.0) then
c following loop added for pc
        do 99 j=1,jr2
   99   call grow('r',11,jr2,g2,nc2,ierr)
        do 100 j=1, -jr1
        jlap=j
        jstart=jstart+1
        jr1=jr1+1
        jr2=jr2+1
cu        call growf4('r', 11, jr2, g2, nc2, ierr)
        call grow('r', 11, jr2, g2, nc2, ierr)
        kc2=ics2-1

        do 50 ilap=1, ncout
        kc2=kc2+1
        if(kc2.gt.ncmax) stop 'kc2 exceeded ncmax'
        mlap1=0
        mlap2=0
        g1lap(ilap,1)=dval
        if(kc2.gt.0) then
          if(g2(kc2).lt.dvtest) then
             mlap1=2
             mlap2=1
             g2lap(ilap,1)=g2(kc2)
             minc=amin0(ilap,minc)
             maxc=amax0(ilap,maxc)
             minr=amin0(jlap,minr)
             maxr=amax0(jlap,maxr)
          else
             g2lap(ilap,1)=dval
          endif
        else
           g2lap(ilap,1)=dval
        endif
        gmap(ilap,1)= float(mlap1)
        gmap(ilap,2)= float(mlap2)
50     continue
cu        call growf4('w', 12, jlap, g1lap(1,1), ncout, ierr)
cu        call growf4('w', 13, jlap, g2lap(1,1), ncout, ierr)
cu        call growf4('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 12, jlap, g1lap(1,1), ncout, ierr)
        call grow('w', 13, jlap, g2lap(1,1), ncout, ierr)
        call grow('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 19, jlap, gmap(1,2), ncout, ierr)
        call grow('w', 16, jlap, g1lap(1,1), ncout, ierr)
100     continue
c       type*,'after 100:jr1,jr2 = ',jr1,jr2
        go to 210
      endif
c 
c case where starting row is outside grid2 only
      if(jrs2.le.0) then
        do 199 j=1,jr1
  199   read(10) dlt
        do 200 j=1, -jr2
        jlap=j
        jstart=jstart+1
        jr1=jr1+1
        jr2=jr2+1
cu        call growf4('r', 10, jr1, g1, nc, ierr)
        call grow('r', 10, jr1, g1, nc, ierr)

        kc1=ics1-1

        do 150 ilap=1, ncout
        kc1=kc1+1
        if(kc1.gt.ncmax) stop 'kc1 exceeded ncmax'
        mlap1=0
        mlap2=0
        g2lap(ilap,1)=dval
        if (nrange(kc1,1).ne.0) then
           g1lap(ilap,1)=dval
        else
           if(g1(kc1).lt.dvtest) then
              mlap1=1
              mlap2=2
              g1lap(ilap,1)=g1(kc1)
              minc=amin0(ilap,minc)
              maxc=amax0(ilap,maxc)
              minr=amin0(jlap,minr)
              maxr=amax0(jlap,maxr)
           else
              g1lap(ilap,1)=dval
          endif
        endif
        gmap(ilap,1)= float(mlap1)
        gmap(ilap,2)= float(mlap2)
150    continue
cu        call growf4('w', 12, jlap, g1lap(1,1), ncout, ierr)
cu        call growf4('w', 13, jlap, g2lap(1,1), ncout, ierr)
cu        call growf4('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 12, jlap, g1lap(1,1), ncout, ierr)
        call grow('w', 13, jlap, g2lap(1,1), ncout, ierr)
        call grow('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 19, jlap, gmap(1,2), ncout, ierr)
        call grow('w', 16, jlap, g2lap(1,1), ncout, ierr)
200     continue
c       type*,'after 200:jr1,jr2 = ',jr1,jr2
      endif
c
c set up for main do-loop
c
210   jstart=jstart+1
      jend=nrout
      if (jre1.gt.nr) jend=nrout-jre1+nr
      if (jre2.gt.nr2) jend=nrout-jre2+nr2
c     type*,'jstart,jend=',jstart,jend
c start main part of do-loop 
      do 400 j=jstart, jend
      jlap=j
      jr1=jr1+1
      jr2=jr2+1
cu      call growf4('r', 10, jr1, g1, nc, ierr)
cu      call growf4('r', 11, jr2, g2, nc2, ierr)
      call grow('r', 10, jr1, g1, nc, ierr)
      call grow('r', 11, jr2, g2, nc2, ierr)
      kc1=ics1-1
      kc2=ics2-1

c     type*
c     type*,'jr1,jr2,jlap=',jr1,jr2,jlap
      do 300 ilap=1, ncout
      kc1=kc1+1
      kc2=kc2+1

      mlap1=0
      mlap2=0
      if(nrange(kc1,1).ne.0) then
           g1lap(ilap,1)=dval
      else
        if(g1(kc1).lt.dvtest) then
           mlap1=1
           g1lap(ilap,1)=g1(kc1)
        else
           g1lap(ilap,1)=dval
        endif
      endif
      if(nrange(kc2,3).ne.0) then
           g2lap(ilap,1)=dval
      else
        if(g2(kc2).lt.dvtest) then
           mlap2=2
           g2lap(ilap,1)=g2(kc2)
        else
           g2lap(ilap,1)=dval
        endif
      endif
      if(g1lap(ilap,1).gt.dvtest.or.g2lap(ilap,1).gt.dvtest) then
        g1lap(ilap,2)=dval
      else
        g1lap(ilap,2)=g1lap(ilap,1)-g2lap(ilap,1)
        ndif=ndif+1
        write(17) dum,dum,g1lap(ilap,2)
      endif
      mlap= (1-mlap1)*mlap2 + mlap1
      mlap2 = (2-mlap2)*mlap1 + mlap2/2
c     type*,'mlap1,mlap2,mlap=',mlap1,mlap2,mlap
      gmap(ilap,1)= float(mlap)
      gmap(ilap,2)= float(mlap2)
      if(mlap.gt.0) then
          minc=amin0(ilap,minc)
          maxc=amax0(ilap,maxc)
          minr=amin0(jlap,minr)
          maxr=amax0(jlap,maxr)
      endif
c
c
c     type*,'ilap,g1lap,g2lap,gmap=',ilap,g1lap(ilap,1),
c    1 g2lap(ilap,1),gmap(ilap,1)
300   continue
cu      call growf4('w', 12, jlap, g1lap(1,1), ncout, ierr)
cu      call growf4('w', 13, jlap, g2lap(1,1), ncout, ierr)
cu      call growf4('w', 14, jlap, gmap(1,1), ncout, ierr)
      call grow('w', 12, jlap, g1lap(1,1), ncout, ierr)
      call grow('w', 13, jlap, g2lap(1,1), ncout, ierr)
      call grow('w', 14, jlap, gmap(1,1), ncout, ierr)
      call grow('w', 19, jlap, gmap(1,2), ncout, ierr)
      call grow('w', 16, jlap, g1lap(1,2), ncout, ierr)
400   continue
c
c write rest of grids if out of range of one of them
c
      if(jend.eq.nrout) go to 666
c case if going outside of grid1 but not grid2
      if(jre1.gt.nr) then
        do 500 j=jend+1, nrout
        jlap=j
        jr1=jr1+1
        jr2=jr2+1
cu        call growf4('r', 11, jr2, g2, nc2, ierr)
        call grow('r', 11, jr2, g2, nc2, ierr)
        kc2=ics2-1

        do 450 ilap=1, ncout
        kc2=kc2+1
        mlap1=0
        mlap2=0
        g1lap(ilap,1)=dval
         if(nrange(kc2,3).ne.0) then
              g2lap(ilap,1)=dval
         else
           if(g2(kc2).lt.dvtest) then
             mlap1=2
             mlap2=1
             g2lap(ilap,1)=g2(kc2)
             minc=amin0(ilap,minc)
             maxc=amax0(ilap,maxc)
             minr=amin0(jlap,minr)
             maxr=amax0(jlap,maxr)
           else
             g2lap(ilap,1)=dval
           endif
         endif
        gmap(ilap,1)= float(mlap1)
        gmap(ilap,2)= float(mlap2)
450     continue
cu        call growf4('w', 12, jlap, g1lap(1,1), ncout, ierr)
cu        call growf4('w', 13, jlap, g2lap(1,1), ncout, ierr)
cu        call growf4('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 12, jlap, g1lap(1,1), ncout, ierr)
        call grow('w', 13, jlap, g2lap(1,1), ncout, ierr)
        call grow('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 19, jlap, gmap(1,2), ncout, ierr)
        call grow('w', 16, jlap, g1lap(1,1), ncout, ierr)
500     continue
c       type*,'after 500:jr1,jr2 = ',jr1,jr2
        go to 666
      endif
c case if going outside of grid2 but not grid1
      if(jre2.gt.nr2) then
        do 600 j=jend+1, nrout
        jlap=j
        jr1=jr1+1
        jr2=jr2+1
cu        call growf4('r', 10, jr1, g1, nc, ierr)
        call grow('r', 10, jr1, g1, nc, ierr)
        kc1=ics1-1

        do 550 ilap=1, ncout
        kc1=kc1+1
        mlap1=0
        mlap2=0
        g2lap(ilap,1)=dval
        if(g1(kc1).lt.dvtest) then
           mlap1=1
           mlap2=2
           g1lap(ilap,1)=g1(kc1)
           minc=amin0(ilap,minc)
           maxc=amax0(ilap,maxc)
           minr=amin0(jlap,minr)
           maxr=amax0(jlap,maxr)
        else
           g1lap(ilap,1)=dval
        endif
        gmap(ilap,1)= float(mlap1)
        gmap(ilap,2)= float(mlap2)
550     continue
        jstart=jstart+1
cu        call growf4('w', 12, jlap, g1lap(1,1), ncout, ierr)
cu        call growf4('w', 13, jlap, g2lap(1,1), ncout, ierr)
cu        call growf4('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 12, jlap, g1lap(1,1), ncout, ierr)
        call grow('w', 13, jlap, g2lap(1,1), ncout, ierr)
        call grow('w', 14, jlap, gmap(1,1), ncout, ierr)
        call grow('w', 19, jlap, gmap(1,2), ncout, ierr)
        call grow('w', 16, jlap, g2lap(1,1), ncout, ierr)
600     continue
c       type*,'after 600:jr1,jr2 = ',jr1,jr2
      endif
c
c overlap grids have been assembled, close grids
666   call gclose(10,'keep')
      call gclose(11,'keep')
      call gclose(16,'keep')
cu      call gclose(12, 'keep')
cu      call gclose(13, 'keep')
cu      call gclose(14, 'keep')
c

      open(15,file='grid1edg.xyz',form='unformatted',status='unknown')
      open(18,file='grid2edg.xyz',form='unformatted',status='unknown')
      call biharm(12,13,14,15,ns1tot,-1.0)
      if(ns1tot.ne.0) then
        call getstats(15,ns1tot,s1ave,s1median,s1mode,ns1m,s1min,s1max)
      endif
      call biharm(13,12,19,18,ns2tot,1.0)
      if(ns2tot.ne.0) then
        call getstats(18,ns2tot,s2ave,s2median,s2mode,ns2m,s2min,s2max)
      endif
      call getstats(17,ndif,dave,dmedian,dmode,ndm,dmin,dmax)
      print 886, igrid1,igrid2,igrid1,igrid2,s1ave,s2ave,dave,s1median,
     1s2median,dmedian
886   format(' Shift statistics:'
c     1 /,10x,'edge of grid1',10x,'edge of grid 2',10x,
     1 /,7x,'edge of ',a12,3x,'edge of 'a12,3x,
     2 a12,' - ',a12,
     3 /,' Mean:    ',g14.7,2(10x,g14.7),
     4 /,' Median:  ',g14.7,2(10x,g14.7))
      print 887, s1mode(1),s2mode(1),dmode(1)
887   format(' Mode(s): ',g14.7,2(10x,g14.7))
888   format(10x,g14.7,2(10x,g14.7))
8881  format(24x,2(10x,g14.7))
8882  format(10x,g14.7,34x,g14.7)
8883  format(58x,g14.7)
      ind=1
900   ns1m=ns1m-1
      ns2m=ns2m-1
      ndm=ndm-1
      ind=ind+1
      if(ns1m.gt.0.and.ns2m.gt.0.and.ndm.gt.0) then
        print 888, s1mode(ind),s2mode(ind),dmode(ind)
      else if(ns1m.gt.0.and.ns2m.gt.0) then
        print 888, s1mode(ind),s2mode(ind)
      else if(ns1m.gt.0.and.ndm.gt.0) then
        print 8882, s1mode(ind),dmode(ind)
      else if(ns2m.gt.0.and.ndm.gt.0) then
        print 8881, s2mode(ind),dmode(ind)
      else if(ns1m.gt.0) then
        print 888, s1mode(ind)
      else if(ns2m.gt.0) then
        print 8881, s2mode(ind)
      else if(ndm.gt.0) then
        print 8883, dmode(ind)
      else
        go to 901
      endif
      go to 900
901   print 893, ns1tot,ns2tot,ndif,s1min,s2min,dmin,s1max,s2max,dmax,
     1igrid1,igrid2
893   format(' #values: ',i10,2(14x,i10),
     1 /,    ' Min:     ',g14.7,2(10x,g14.7),
     2 /,    ' Max:     ',g14.7,2(10x,g14.7),
     3//' Subtract shift from ',a12,' or add to ',a12
     4//' Resolve multiple modes by examining the histogram of g1mg2lap.
     5grd using program',/,' ghist or of grid1edg.xyz or grid2edg.xyz us
     6ing program phist.')
      call gclose(12,'keep')
      call gclose(13,'keep')
      call gclose(14,'delete')
      close(15)
      close(17)
      close(18)
      call gclose(19,'delete')
c      stop
      end

c**************************************************************************
      subroutine biharm(iunit,junit,lunit,kunit,nstotal,sign)
c iunit contains grid1 in the overlap, trimmed to define the merge boundary
c junit contains grid2 in the overlap area
c kunit is the output xyz file of shifts along the merge boundary
c lunit is the data coverage grid
c nstotal is the total number of shifts computed
c sign is used to adjust the sign of the output shifts
      parameter (ncmax=3000)
      dimension s(ncmax)
c      dimension shft(2*ncmax)
      dimension u(13), w(13)
      character title3*56
      common /data/ g1lap(ncmax,5),g2lap(ncmax,5),gmap(ncmax,5)
      common /bih/ title3,dval,ncout,xoout,yoout,minc,maxc,minr,maxr
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      common /sums/ nsum1,nsum2
c
      if(sign.lt.0.0) then
        sign = -1.0
      else
        sign = 1.0
      endif
      ierr=0
cu      call gclose(12, 'keep')
cu      call gclose(13, 'keep')
cu      call gclose(14, 'keep')
cu      call gopen(12, 'grid1lap.gd', 'old', 'read', ierr)
cu      call gopen(13, 'grid2lap.gd', 'old', 'read', ierr)
cu      call gopen(14, 'maplap.gd', 'old', 'read', ierr)
cu      call gh1f4('r',12,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
cu      call gh1f4('r',13,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
cu      call gh1f4('r',14,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      call gheader('r',iunit,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      call gheader('r',junit,title3,ncout,nrout,xoout,dx,yoout,dy,ierr)
      rewind(lunit)
c
c  initialize s array
      do ii=1,ncout
      s(ii) = dval
      enddo
c
c  as look at shifts, keep track of min, max
      smin=dval
      smax=-dval
      nstotal=0
c
c first time through--read first 5 rows beginning at minr
      krow=minr-1
      do j=1,5
      indx=j
      krow=krow+1
cu      call growf4('r', 12, krow, g1lap(1,indx), ncout, ierr)
cu      call growf4('r', 13, krow, g2lap(1,indx), ncout, ierr)
cu      call growf4('r', 14, krow, gmap(1,indx), ncout, ierr)
      call grow('r', iunit, krow, g1lap(1,indx), ncout, ierr)
      call grow('r', junit, krow, g2lap(1,indx), ncout, ierr)
      call grow('r', lunit, krow, gmap(1,indx), ncout, ierr)
      enddo
      ypos=yoout+2*dx
      go to 700
c
c other times through--shift arrays up one row and read new 5th row
c
650   ypos=ypos+dy
      do jk=1,4
        do ii=minc,maxc
        g1lap(ii,jk)=g1lap(ii,jk+1)
        g2lap(ii,jk)=g2lap(ii,jk+1)
        gmap(ii,jk)=gmap(ii,jk+1)
        enddo
      enddo
cu      call growf4('r', 12, krow, g1lap(1,5), ncout, ierr)
cu      call growf4('r', 13, krow, g2lap(1,5), ncout, ierr)
cu      call growf4('r', 14, krow, gmap(1,5), ncout, ierr)
      call grow('r', iunit, krow, g1lap(1,5), ncout, ierr)
      call grow('r', junit, krow, g2lap(1,5), ncout, ierr)
      call grow('r', lunit, krow, gmap(1,5), ncout, ierr)
c
c compute shift for each col along 3rd row of 5-row arrays
c  check to see if enough data points for computation using
c  nsum1 (no. of pts in grid1) and nsum2 (same for grid2) 
c  if data for u(7) not available for grid1, go to next
c
700   do 750 icol=minc+2,maxc-2
c      print*,gmap(icol,3)
      if(gmap(icol,3).ne.1.0) go to 725
      xpos=xoout+(icol-1)*dx
      nsum1 = 0
      nsum2 = 0
c
c determine u(1) and weight w(1)
      ic=icol
      call usetup(5, ic, u(1), w(1))
c
c determine u(2) through u(4) and u(10) through u(12) and wts
      ic=icol-1
      do mm=2,4
      m=mm
      call usetup(4, ic, u(m), w(m))
      call usetup(2, ic, u(8+m), w(8+m)) 
      ic=ic + 1
      enddo
c
c determine u(5) through u(9) and wts
      ic=icol-2
      do mm=5,9
      m=mm
      call usetup(3, ic, u(m), w(m))
      ic=ic+1
      enddo
c
c determine u(13) and wt(13)
      ic=icol
      call usetup(1, ic, u(13), w(13))
c
c If nsum1 + nsum2 = 13 then all data are valid.
c  if nsum1 & 2 both >=4, then there are at least 4 points in each grid
c  if not, go to next computation
c
c     type*,'krow,icol,nsum1,nsum2=',krow,icol,nsum1,nsum2
      if(nsum1+nsum2.lt.13.or.nsum1.lt.4.or.nsum2.lt.4) go to 725
c 
c evaluate the biharmonic operator and find constant shift
c
c      print*,u,nsum1,nsum2,gmap(icol,3)
      bih=20.0e0*u(7)
      bih=bih - 8.0e0*(u(3) + u(8) + u(11) + u(6))
      bih=bih + 2.0e0*(u(2) + u(4) + u(12) + u(10))
      bih=bih + u(1) + u(9) + u(13) + u(5)
      wtsumw2= -8.0e0*(w(3) + w(8) + w(11) + w(6))
      wtsumw2=wtsumw2 + 2.0e0*(w(2) + w(4) + w(12) + w(10))
      wtsumw2=wtsumw2 + w(1) + w(9) + w(13) + w(5)
      nstotal=nstotal+1
c      if(nstotal.gt.2*ncmax) stop 'shift array too small'
      shft = sign*bih/wtsumw2
      write(kunit) xpos,ypos,shft
      s(icol)=shft
      ssum=ssum+shft
      smin= amin1 (shft,smin)
      smax= amax1 (shft,smax)
      go to 750
725   s(icol)=dval
750   continue      
c
c write out the s array for row=krow-2
c
      ksrow=krow-2
c next row
      krow=krow+1
      if (krow.le.maxr) go to 650
c end of row loop
c
c      if(nstotal.eq.0) then
c        stop 'No points could be checked.  Surveys not touching?'
c      endif
      return
      end
c**************************************************************************
      subroutine usetup (jr, ic, uval, w2val) 
      parameter (ncmax=3000)
      common /data/ g1lap(ncmax,5),g2lap(ncmax,5),gmap(ncmax,5)
      common /sums/ nsum1, nsum2
c this subroutine determines the value of u and wt for grid2
      gmval=gmap(ic,jr)
      ichk=nint( alog( gmval+1.0 ))
      chk=float(ichk)
      gmap1=  2.0 - gmval
      gmap2=  gmval - 1.0
c gmap1=1 and gmap2=0 to use value from grid1
c gmap1=0 and gmap2=1 to use value from grid2
c if no data avail (chk=0), nsum's are made zero.
      uval= gmap1*g1lap(ic,jr) + gmap2*g2lap(ic,jr)
      w2val=gmap2
      nsum1 = (nsum1 + gmap1) * chk
      nsum2 = (nsum2 + gmap2) * chk
      return
      end

c***************************************************************************
      subroutine rcfind(xyo, del, f, nf)

      ff= (f - xyo)/del + 1.e0
      nf= nint( ff )

      return
      end
c****************************************************************************
c function nrange checks to see if a given row or col is within the
c range of the grid specs
      function nrange(n,icode)
c icode
c   1     check of col range in grid1
c   2     check of row range in grid1
c   3     check of col range in grid2
c   4     check of row range in grid2
      common /specs/xo,yo,nc,nr,dx,dy,xo2,yo2,nc2,nr2
      nrange=0
      if (n. le. 0) then
         nrange=1
         return
      endif
      if (icode.eq.1) ncheck=nc
      if (icode.eq.2) ncheck=nr
      if (icode.eq.3) ncheck=nc2
      if(icode.eq.4) ncheck=nr2
      if(n .gt. ncheck) nrange=1
      return
      end
c****************************************************************************
c function amode works like amod should
      function amode(x,y)
      if(y.eq.0.0) then
        amode=0.0
      else
        amode=x-anint(x/y)*y
      endif
      return
      end
c****************************************************************************
      subroutine pmode(n, nm, dm, dmed)
      parameter (ncmax=3000)
      dimension dm(100)
      common /data/a(ncmax*15)
      nm=1
c bubble sort a() in ascending order
c      if(n.lt.2) return
c      do 21 i=1,n-1
c      do 21 j=i+1,n
c      if(a(j)-a(i)) 20,21,21
c   20 ax=a(i)
c      a(i)=a(j)
c      a(j)=ax
c   21 continue
      call quicksort(n)
c get median
      if(2*(n/2).eq.n) then
        dmed=0.5*(a(n/2)+a(n/2+1))
      else
        dmed=a((n+1)/2)
      endif
c get mode(s)
      aminold=-1.e38
      amaxold=1.e38
c set initial window size to n/2, decimate at each step
      j=n
   29 j=j/2
      if(j.lt.3) return
c locate the maximum of the distribution for this window size
      pmax=-1.e38
      imax=0
      do 30 i=1,n-j
      den=n*(a(i+j)-a(i))
      if(den.eq.0.) go to 30
      p=float(j)/den
      if(p.gt.pmax) then
        pmax=p
        imax=i
      endif
   30 continue
      amin=a(imax)
      amax=a(imax+j)
      dmode=0.5*(amin+amax)
c is this maximum consistent with results from the previous window size?
      if(dmode.lt.aminold.or.dmode.gt.amaxold) then
c no, we have a new mode
c increment nm and save the new mode
         nm=nm+1
         if(nm.gt.100) stop 'too many modes'
         dm(nm)=dmode
      else
c yes, update the old mode with the new value
        dm(nm)=dmode
      endif
      aminold=amin
      amaxold=amax
      go to 29
      end
c****************************************************************************
      subroutine quicksort(n)
c implementation of the quicksort algorithm
      parameter (ncmax=3000)
      integer hi(11),lo(11)
      common /data/a(ncmax*15)
      lo(1)=1
      hi(1)=n
      is=1
      do while(is.gt.0)
        if(lo(is).ge.hi(is)) then
          is=is-1
          cycle
        endif
        i=lo(is)
        j=hi(is)
        pivot=a(j)
        do while(i.lt.j)
          do while((i.lt.j).and.(a(i).le.pivot))
            i=i+1
          enddo
          do while((j.gt.i).and.(a(j).ge.pivot))
            j=j-1
          enddo
          if(i.lt.j) then
            temp=a(i)
            a(i)=a(j)
            a(j)=temp
          endif
        enddo
        j=hi(is)
        temp=a(i)
        a(i)=a(j)
        a(j)=temp
        if((i-lo(is)).lt.(hi(is)-i)) then
          lo(is+1)=lo(is)
          hi(is+1)=i-1
          lo(is)=i+1
        else
          lo(is+1)=i+1
          hi(is+1)=hi(is)
          hi(is)=i-1
        endif
        is=is+1
      enddo
      return
      end
c****************************************************************************
      subroutine getstats(iunit,nd,dave,dmedian,dmode,nmode,dmin,dmax)
      parameter (ncmax=3000)
      real*8 ave8
      dimension dmode(2,100)
      common /data/glap(ncmax*15)
      if(nd.gt.ncmax*15) then
        print*,'nd = ',nd,' > ncmax*15 = ',ncmax*15
        stop 'Cannot compute statistics'
      endif
      rewind(iunit)
      dmin=1.e38
      dmax=-1.e38
      ave8=0.
c read the data off iunit, compute the min, max, and mean
      do 10 i=1,nd
      read(iunit) x,y,glap(i)
      ave8=ave8+glap(i)
      if(glap(i).lt.dmin) dmin=glap(i)
      if(glap(i).gt.dmax) dmax=glap(i)
   10 continue
      dave = ave8/dfloat(nd)
c compute the median amd mode(s)
      call pmode(nd,nmode,dmode,dmedian)
      return
      end
