program gsplat

*** overwrite a "master" grid file with a second grid file
*** grids need not conform in boundary or granularity (kinds must match)
*** output grid boundary/granularity identical to "master"

      implicit double precision(a-h,o-z)

c - hack!!!! changed len to 5001*5001

      parameter(len=14500*20500)
c     parameter(len=3001*3841)
c      parameter(len=2049*2049)
      character*63 fname
      logical makpos,flag1
      integer h(len),hrec(18001)
      real*4  z(len),zrec(18001)
      equivalence(h(1),z(1)),(hrec(1),zrec(1))

      lin=1
      lout=2

*** load the master grid

      write(6,1)
    1 format('program gsplat    --> Enter master file name: ',$)
      read(5,'(a)') fname
      open(lin,file=fname,status='old',form='unformatted')
      read(lin) glamn,glomn,dgla,dglo,nla,nlo,ikind
      glamx=glamn+dgla*(nla-1)

      flag1=makpos(glomn)
      if(flag1.and.dglo.lt.0.d0) dglo=-dglo
      glomx=glomn+dglo*(nlo-1)

      if(nla*nlo.gt.len) stop 12345
      if(ikind.eq.0) then
        call r1(lin,h,nla,nlo)
      else
        call r2(lin,z,nla,nlo)
      endif
      close(lin)

*** now read the secondary grid
*** it will overwrite master grid values whenever possible

      write(6,2)
    2 format('                  --> Enter second file name: ',$)
      read(5,'(a)') fname
      open(lin,file=fname,status='old',form='unformatted')
      read(lin) glamn2,glomn2,dgla2,dglo2,nla2,nlo2,ikind2
      if(ikind.ne.ikind2) then
        write(*,*) 'kinds do not match'
        stop 66666
      endif

      flag1=makpos(glomn2)
      if(flag1.and.dglo2.lt.0.d0) dglo2=-dglo2

      icount=0
      iclip=0
      if(ikind.eq.0) then
        do 10 i=1,nla2
          gla=glamn2+(i-1)*dgla2
          read(lin) (hrec(j),j=1,nlo2)
          do 11 j=1,nlo2
            glo=glomn2+(j-1)*dglo2
            icount=icount+1
            call put1(hrec(j),gla,glo,
     *                h,nla,nlo,glamn,dgla,glomn,dglo,iclip)
   11     continue
   10   continue
      else
        do 20 i=1,nla2
          gla=glamn2+(i-1)*dgla2
          read(lin) (zrec(j),j=1,nlo2)
          do 21 j=1,nlo2
            glo=glomn2+(j-1)*dglo2
            icount=icount+1
            call put2(zrec(j),gla,glo,
     *                z,nla,nlo,glamn,dgla,glomn,dglo,iclip)
   21     continue
   20   continue
      endif

*** write the updated grid

      write(6,3)
    3 format('      grid output --> Enter file name: ',$)
      read(5,'(a)') fname
      open(lout,file=fname,status='new',form='unformatted')

      write(lout) glamn,glomn,dgla,dglo,nla,nlo,ikind
      if(ikind.eq.0) then
        call w1(lout,h,nla,nlo)
      else
        call w2(lout,z,nla,nlo)
      endif

      write(6,*) ' icount, iclip = ',icount,iclip

      stop
      end
      subroutine put1(ival,gla,glo,h,nla,nlo,glamn,dgla,glomn,dglo,iclp)

*** load value into grid

      implicit double precision(a-h,o-z)
      integer h(nla,nlo)

      i=idnint((gla-glamn)/dgla)+1
      if(i.lt.1.or.i.gt.nla) then
        iclp=iclp+1
        return
      endif
      j=idnint((glo-glomn)/dglo)+1
      if(j.lt.1.or.j.gt.nlo) then
        iclp=iclp+1
        return
      endif
      h(i,j)=ival

      return
      end
      subroutine put2(val,gla,glo,z,nla,nlo,glamn,dgla,glomn,dglo,iclip)

*** load value into grid

      implicit double precision(a-h,o-z)
      real*4 z(nla,nlo),val

      i=idnint((gla-glamn)/dgla)+1
      if(i.lt.1.or.i.gt.nla) then
        iclip=iclip+1
        return
      endif
      j=idnint((glo-glomn)/dglo)+1
      if(j.lt.1.or.j.gt.nlo) then
        iclip=iclip+1
        return
      endif
      z(i,j)=val

      return
      end
      subroutine w1(lout,h,nla,nlo)

*** write records south to north (elements are west to east)

      implicit double precision(a-h,o-z)
      integer h(nla,nlo)

      do 1 i=1,nla
    1 write(lout) (h(i,j),j=1,nlo)

      return
      end
      subroutine w2(lout,z,nla,nlo)

*** write records south to north (elements are west to east)

      implicit double precision(a-h,o-z)
      real*4 z(nla,nlo)

      do 1 i=1,nla
    1 write(lout) (z(i,j),j=1,nlo)

      return
      end
      logical function makpos(glon)

*** insure longitude is positive

      implicit double precision(a-h,o-z)

      makpos=.false.

    1 if(glon.lt.0.d0) then
        glon=glon+360.d0
        makpos=.true.
        go to 1
      endif

    2 if(glon.ge.360.d0) then
        glon=glon-360.d0
        go to 2
      endif

      return
      end
      subroutine r1(lin,h,nla,nlo)

***  read records south to north (elements are west to east)

      implicit double precision(a-h,o-z)
      integer h(nla,nlo)

      do 1 i=1,nla
    1 read(lin) (h(i,j),j=1,nlo)

      return
      end
      subroutine r2(lin,z,nla,nlo)

***  read records south to north (elements are west to east)

      implicit double precision(a-h,o-z)
      real*4 z(nla,nlo)

      do 1 i=1,nla
    1 read(lin) (z(i,j),j=1,nlo)

      return
      end