c--------( GMERGE : merge multiple grid data )--------
      parameter (LX=4051, LY=4051)
c
      dimension f(LX,LY), g(LX,LY), r(LX,LY)
      dimension a(LX,LY), h(LX,LY), t(LX,LY)
      character*80 wdr, flog, fnam1, fnam2
      character buf*80, area*8, narea*8, area1*8, anul*10, salt*8
c
      call premsg('GMERGE : merge multiple grid data')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnam1 = wdr
      fnam2 = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call premsg('')
      call gparma('Enter Name-label for new grid ==> ', 8, narea)
      call gparmi('Select Projection-Number ==> ', nccn)
      if(nccn.ge.200) ncn = nccn - 200
      iorgn = 35*60
      korgn = 135*60
      ispan = 0
      ispbn = 0
      if(ncn.ge.1.and.ncn.le.62) then
        iorgn = 0
        korgn = (ncn-30)*360 - 180
        if(ncn.gt.60) korgn = 0
        if(ncn.eq.61) iorgn = +5400
        if(ncn.eq.62) iorgn = -5400
      else if(ncn.eq.72) then
        ispan = 30*60
        ispbn = 40*60
      endif
      call premsg('Specify gridding parameters')
      call premsg('  Southwest corner Coord. in km ?')
      call gparmf('    Northing : ', xl)
      ixl = nint(xl*1000.)
      call gparmf('    Easting  : ', yl)
      iyl = nint(yl*1000.)
      call gparmi('  Unit Mesh size in m :  ', msz)
      if(msz.le.0) call abendm('invalid value')
      call gparmi('    Number of Mesh to N :  ', mxn)
      if(mxn.gt.LX) call abendm('too large array size')
      call gparmi('    Number of Mesh to E :  ', myn)
      if(myn.gt.LY) call abendm('too large array size')
      vnuln = 9999.9
      call gparmf('Boundary zone Width ?  (in km) ==> ', sb)
      nb = nint(sb*1000./float(msz))
      mxnb = mxn + nb + nb
      mynb = myn + nb + nb
      if(mxnb.gt.LX .or. mynb.gt.LY)
     *     call abendm('too large array size')
      dnorm = float(nb)
      nddx = nb*nb*4
      call gparma('Enter new grid output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      open(2,file=fnam2,status='new')
      do 100 j=1,mynb
      do 100 i=1,mxnb
        g(i,j) = vnuln
        h(i,j) = vnuln
  100 continue
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'GMERGE : merge multiple grid data'       
      write(6,'(a,a)') ' Output new file : ', fnam2
      write(6,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km'
      write(6,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *              xl, yl, msz, msz, mxn, myn, '    var.'
      open(9,status='scratch')
c
      nsrc = 0
    1 nsrc = nsrc + 1
      call premsg('')
      call premsg('Enter Input source data filename')
      call gparma('             (Hit RETURN to End) ==> ',
     *            81-ldr, fnam1(ldr:80))
      if(fnam1(ldr:ldr).eq.' ') goto 8
      open(1,file=fnam1,status='old')
    2 read(1,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      endif
      if(ncc.ne.nccn .or.
     *   iorg.ne.iorgn.or.korg.ne.korgn .or.
     *   ispa.ne.ispan.or.ispb.ne.ispbn)
     *     call abendm('Header inconsistent')
      read(1,'(a)') buf
      read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt
      write(salt,'(1x,f7.0)') alt
      if (alt.eq.0.) salt = '    var.'
      if (alt.lt.0.) salt = '   undef'
      if (mszx.ne.msz.or.mszy.ne.msz) call abendm('illegal mesh-size')
      is = (ixs-ixl)/mszx
      js = (iys-iyl)/mszy
      if (is*mszx+ixl.ne.ixs .or.
     *    js*mszy+iyl.ne.iys) call abendm('ranges conflict')
      if(mx.gt.LX.or.my.gt.LY) call abendm('too large array size')
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
      read(1,*) ((r(i,j),i=1,mx),j=1,my)
      if(alt.eq.0.) goto 3
      if(alt.lt.0.) then
        call gparmf('  Altitude [in m] for this data ==> ', alt)
        if (alt.le.0.) call abendm('invalid value')
        write(salt,'(1x,f7.0)') alt
      endif
      goto 4
    3 read(1,'(a)',end=4) buf
      if(buf(1:1).eq.'#') goto 3
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      nc1 = ncc1
      if(ncc1.ge.200) nc1 = ncc1 - 200
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if(area1.ne.area.or.ncc1.ne.ncc .or.
     *   iorg1.ne.iorg.or.korg1.ne.korg .or.
     *   ispa1.ne.ispa.or.ispb1.ne.ispb)
     *     call abendm('Header1 inconsistent')
      read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
      if(ixs1.ne.ixs.or.iys1.ne.iys .or.
     *   mszx1.ne.mszx.or.mszy1.ne.mszy .or.
     *   mx1.ne.mx.or.my1.ne.my)
     *     call abendm('Header2 inconsistent')
      write(anul,'(f10.3)') vnul1
      if(anul(9:10).ne.'00'.or.anul(1:1).ne.' ')
     *     call abendm('unable to handle vnul1 value')
      read(1,*) ((t(i,j),i=1,mx),j=1,my)
    4 close(1)
      write(6,'(a,a,a8,a,i4,4x,4i8)') '--- ',
     *     'Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      write(6,'(4x,a,2a10,2x,2a8,2a6,a7)')
     *     '  (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *     'alt'
      write(6,'(4x,a,2f10.2,2i8,2i6),a8')
     *     'source data ', xs, ys, mszx, mszy, mx, my, salt
      write(9,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb
      write(9,*) xs, ys, mszx, mszy, mx, my, alt
      do 200 j=1,my
        jp = js + j + nb
        if (jp.le.0.or.jp.gt.mynb) goto 200
        do 210 i=1,mx
          if (r(i,j).eq.vnul) goto 210
          ip = is + i + nb
          if (ip.le.0.or.ip.gt.mxnb) goto 210
          g(ip,jp) = r(i,j)
          if (alt.ne.0.) then
            h(ip,jp) = alt
          else if (t(i,j).ne.vnul1) then
            h(ip,jp) = t(i,j)
          else
            write(buf,'(a30,f8.1,a8,i4,a4,i4)')
     *             '  t(i,j)=vnul  where  r(i,j) =', r(i,j),
     *             '  at  i=', i, '  j=', j
            call premsg(buf(1:58))
            call abendm('Aborted')
          endif
  210 continue
  200 continue
      if (nsrc.eq.1) then
        do 300 j=1,mynb
        do 300 i=1,mxnb
          f(i,j) = g(i,j)
          a(i,j) = h(i,j)
          g(i,j) = vnuln
          h(i,j) = vnuln
  300   continue
      else
        call prompt('... Gmerge processing : ')
        call dpcent(0, 0)
        do 400 j=1,mynb
        lf = 0
        lg = 0
        do 400 i=1,mxnb
          ndd = nddx
          if (lf.lt.0) then
            do 500 jm=-nb,+nb
              if ((i+nb).gt.mxnb) goto 500
              if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 500
              if (f(i+nb,j+jm).ne.vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp.lt.ndd) ndd = nddp
              endif
  500       continue
            if (ndd.ne.nddx) then
              lf = 0
              df = -sqrt(float(ndd)) / dnorm
            endif
          else if (lf.gt.0) then
            do 510 jm=-nb,+nb
              if ((i+nb).gt.mxnb) goto 510
              if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 510
              if (f(i+nb,j+jm).eq.vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp.lt.ndd) ndd = nddp
              endif
  510       continue
            if (ndd.ne.nddx) then
              lf = 0
              df = sqrt(float(ndd-1)) / dnorm
            endif
          else
            if (f(i,j).eq.vnuln) then
              do 600 jm=-nb,+nb
              do 600 im=-nb,+nb
                if ((i+im).le.0.or.(i+im).gt.mxnb) goto 600
                if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 600
                if (f(i+im,j+jm).ne.vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp.lt.ndd) ndd = nddp
                endif
  600         continue
              if (ndd.eq.nddx) then
                lf = -1
                df = -2.
              else
                df = -sqrt(float(ndd)) / dnorm
              endif
            else
              do 610 jm=-nb,+nb
              do 610 im=-nb,+nb
                if ((i+im).le.0.or.(i+im).gt.mxnb) goto 610
                if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 610
                if (f(i+im,j+jm).eq.vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp.lt.ndd) ndd = nddp
                endif
  610         continue
              if (ndd.eq.nddx) then
                lf = +1
                df = +2.
              else
                df = sqrt(float(ndd-1)) / dnorm
              endif
            endif
          endif
          ndd = nddx
          if (lg.lt.0) then
            do 700 jm=-nb,+nb
              if ((i+nb).gt.mxnb) goto 700
              if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 700
              if (g(i+nb,j+jm).ne.vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp.lt.ndd) ndd = nddp
              endif
  700       continue
            if (ndd.ne.nddx) then
              lg = 0
              dg = -sqrt(float(ndd)) / dnorm
            endif
          else if (lg.gt.0) then
            do 710 jm=-nb,+nb
              if ((i+nb).gt.mxnb) goto 710
              if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 710
              if (g(i+nb,j+jm).eq.vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp.lt.ndd) ndd = nddp
              endif
  710       continue
            if (ndd.ne.nddx) then
              lg = 0
              dg = sqrt(float(ndd-1)) / dnorm
            endif
          else
            if (g(i,j).eq.vnuln) then
              do 800 jm=-nb,+nb
              do 800 im=-nb,+nb
                if ((i+im).le.0.or.(i+im).gt.mxnb) goto 800
                if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 800
                if (g(i+im,j+jm).ne.vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp.lt.ndd) ndd = nddp
                endif
  800         continue
              if (ndd.eq.nddx) then
                lg = -1
                dg = -2.
              else
                dg = -sqrt(float(ndd)) / dnorm
              endif
            else
              do 810 jm=-nb,+nb
              do 810 im=-nb,+nb
                if ((i+im).le.0.or.(i+im).gt.mxnb) goto 810
                if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 810
                if (g(i+im,j+jm).eq.vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp.lt.ndd) ndd = nddp
                endif
  810         continue
              if (ndd.eq.nddx) then
                lg = +1
                dg = +2.
              else
                dg = sqrt(float(ndd-1)) / dnorm
              endif
            endif
          endif
          if (df.lt.0.) then
            if (dg.lt.0.) then
              if ((1.+df+dg).gt.0.) then
                avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm)
                avg = avrg(g,i,j,mxnb,mynb,vnuln, (1.+df-dg)*dnorm)
                r(i,j) = (avf*(1.+df) + avg*(1.+dg)) / (2.+df+dg)
                ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm)
                avh = avrg(h,i,j,mxnb,mynb,vnuln, (1.+df-dg)*dnorm)
                t(i,j) = (ava*(1.+df) + avh*(1.+dg)) / (2.+df+dg)
              else
                r(i,j) = vnuln
                t(i,j) = vnuln
              endif
            else if (dg.lt.1.) then
              if ((1.+df-dg).gt.0.) then
                avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
                avg = avrg(g,i,j,mxnb,mynb,vnuln, (1.+df-dg)*dnorm)
                r(i,j) = (avf*(1.+df-dg) + avg*(1.+dg)) / (2.+df)
                ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
                avh = avrg(h,i,j,mxnb,mynb,vnuln, (1.+df-dg)*dnorm)
                t(i,j) = (ava*(1.+df-dg) + avh*(1.+dg)) / (2.+df)
              else
                r(i,j) = g(i,j)
                t(i,j) = h(i,j)
              endif
            else
              r(i,j) = g(i,j)
              t(i,j) = h(i,j)
            endif
          else if (df.lt.1.) then
            if (dg.lt.0.) then
              if ((1.+dg-df).gt.0.) then
                avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm)
                avg = avrg(g,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
                r(i,j) = (avf*(1.+df) + avg*(1.+dg-df)) / (2.+dg)
                ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm)
                avh = avrg(h,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
                t(i,j) = (ava*(1.+df) + avh*(1.+dg-df)) / (2.+dg)
              else
                r(i,j) = f(i,j)
                t(i,j) = a(i,j)
              endif
            else if (dg.lt.1.) then
              avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
              avg = avrg(g,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
              r(i,j) = (avf*(1.+df-dg) + avg*(1.+dg-df)) / 2.
              ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
              avh = avrg(h,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
              t(i,j) = (ava*(1.+df-dg) + avh*(1.+dg-df)) / 2.
            else
              avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
              r(i,j) = (avf*df + g(i,j)*(2.-df)) / 2.
              ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm)
              t(i,j) = (ava*df + h(i,j)*(2.-df)) / 2.
            endif
          else
            if (dg.lt.0.) then
              r(i,j) = f(i,j)
              t(i,j) = a(i,j)
            else if (dg.lt.1.) then
              avg = avrg(g,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
              r(i,j) = (f(i,j)*(2.-dg) + avg*dg) / 2.
              avh = avrg(h,i,j,mxnb,mynb,vnuln, (1.-dg)*dnorm)
              t(i,j) = (a(i,j)*(2.-dg) + avh*dg) / 2.
            else
              r(i,j) = (f(i,j) + g(i,j)) / 2.
              t(i,j) = (a(i,j) + h(i,j)) / 2.
            endif
          endif
          call dpcent((j-1)*mxnb+i, mxnb*mynb)
  400   continue
        do 900 j=1,mynb
        do 900 i=1,mxnb
          f(i,j) = r(i,j)
          a(i,j) = t(i,j)
          g(i,j) = vnuln
          h(i,j) = vnuln
  900   continue
      endif
      goto 1
c
    8 call clspin()
      endfile(9)
      rewind(9)
      write(2,'(a)') '# prog.GMERGE output'
      write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
      write(2,'(2i12,4i6,1x,f7.1, 1x,f7.0)')
     *     ixl,iyl, mszx,mszy, mxn,myn, vnuln, 0.
      do 50 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (f(i+nb,j+nb),i=1,mxn)
   50 continue
      write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
      write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
     *     ixl,iyl, mszx,mszy, mxn,myn, vnuln, -1.
      do 60 j=1,myn
        write(2,'((f7.1,9(1x,f7.1)))') (a(i+nb,j+nb),i=1,mxn)
   60 continue
      close(2)
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        write(99,'(//a)') 'GMERGE : merge multiple grid data'
        write(99,'(a,a)') ' Output new file : ', fnam2
        write(99,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km'
        write(99,1000) narea, nccn, iorgn,korgn, ispan,ispbn,
     *                xl, yl, msz, msz, mxn, myn, '    var.'
   92   read(9,'(a8,i4,4x,4i8)',end=93)
     *       area, ncc, iorg, korg, ispa, ispb
        read(9,*) xs, ys, mszx, mszy, mx, my, alt
        write(99,'(a,a,a8,a,i4,4x,4i8)') '--- ',
     *       'Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
        write(99,'(4x,a,2a10,2x,2a8,2a6)')
     *       '  (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my',
     *       'alt'
        write(99,'(4x,a,2f10.2,2i8,2i6,a8)')
     *       'source data ', xs, ys, mszx, mszy, mx, my, salt
        goto 92
   93   close(99)
      endif
      close(9)
c
      stop
 1000 format('Areaname : ', a8, 3x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'gen.new grid', 2f10.2, 2i8, 2i6, a8)
      end
c
      function avrg(f,i,j,mx,my,vnul, d)
      parameter (LX=4051, LY=4051)
      dimension f(LX,LY)
      dd = d * d
      n = nint(d)
      s = 0.
      k = 0
      do 10 jm=-n,+n
        if((j+jm).le.0.or.(j+jm).gt.my) goto 10
        do 20 im=-n,+n
          if((i+im).le.0.or.(i+im).gt.mx) goto 20
          if (f(i+im,j+jm).eq.vnul) goto 20
          ddp = float(im*im + jm*jm)
          if (ddp.gt.dd) goto 20
          s = s + f(i+im,j+jm)
          k = k + 1
   20 continue
   10 continue
      if (k.eq.0) call abendm('Error in function subprogram (argv)')
      avrg = s / float(k)
      return
      end
