!--------( GMERGE : merge multiple grid data )--------
  implicit none
  integer :: lwkdir, lrtrim;  real :: avrg
  real,allocatable :: f(:,:), g(:,:), r(:,:), a(:,:), h(:,:), t(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, buf
  character(20) :: sdtm;  character(10) :: anul
  character(8) :: area, narea, area1, salt
  integer :: ldr, i, j, i1, i2, nb, mxnb, mynb, nddx, ndd, nddp, IER
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ncn, nccn, iorgn, korgn, ispan, ispbn
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixl, iyl, msz, mxn, myn
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: nsrc, is, js, ip, jp, lf, lg, im, jm
  real :: vnul, alt, xs, ys, vnuln, sb, dnorm, xl, yl, vnul1
  real :: df, dg, avf, avg, ava, avh

  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 Coordinate Number (+800 for Old-TOKYO) ==> ', nccn)
  ncn = nccn
  if (nccn >= 800) ncn = nccn - 800
  ispan = 0;  ispbn = 0
  if ((ncn >= 1) .and. (ncn <= 62)) then
    iorgn = 0
    korgn = (ncn-30)*360 - 180
    if (ncn > 60) korgn = 0
    if (ncn == 61) iorgn = +5400
    if (ncn == 62) iorgn = -5400
  else if (ncn == 65) then
    iorgn = 0
    call gparmi2('  Central Longitude ?  (Deg. Min.) : ', i1,i2)
    korgn = i1*60 + i2
  else
    if (ncn == 72) call abendm('nc=72 not implemented')
    call gparmi2('  Origin Latitude ?   (Deg. Min.) : ', i1,i2)
    iorgn = i1*60 + i2
    call gparmi2('  Origin Longitude ?  (Deg. Min.) : ', i1,i2)
    korgn = i1*60 + i2
  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 <= 0) call abendm('invalid value')
  call gparmi('    Number of Mesh to N :  ', mxn)
  call gparmi('    Number of Mesh to E :  ', myn)
  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

  allocate(f(mxnb,mynb));  allocate(g(mxnb,mynb));  allocate(r(mxnb,mynb))
  allocate(a(mxnb,mynb));  allocate(h(mxnb,mynb));  allocate(t(mxnb,mynb))

  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')
  g(:,:) = vnuln;  h(:,:) = vnuln
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'GMERGE : merge multiple grid data'       
  write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km'
  write(6,'(a12,a8,6x,a6,i4,4x,4i8/)')   &
    'Areaname :  ', narea, 'Coord.', nccn, iorgn, korgn, ispan, ispbn
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'gen.new grid', xl, yl, msz, msz, mxn, myn, '    var.'
  call strdtm(sdtm)
  open(9,status='scratch')

  nsrc = 0
  do
    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) == ' ') exit

    open(1,file=fnam1,status='old')
    read(1,'(a)') buf
    do while (buf(1:1) == '#')
      read(1,'(a)') buf
    enddo
    if (buf(13:16) == '    ') then      !-- v2005
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
      if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
      nc = ncc
      if (ncc < 200) then
        nc = ncc; ncc = nc + 800
      else
        ncc = ncc - 200; nc = ncc
      endif
    else                                !-- v2018
      read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
      nc = ncc
      if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
    endif
    if ((nc >= 1) .and. (nc <= 62)) then
      iorg = 0;  ispa = 0;  ispb = 0
      korg = (nc-30)*360 - 180
      if (nc > 60) korg = 0
      if (nc == 61) iorg = +5400
      if (nc == 62) iorg = -5400
    endif
    if ((ncc /= nccn) .or. (iorg /= iorgn) .or. (korg /= korgn)   &
        .or. (ispa /= ispan) .or. (ispb /= 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 == 0.) salt = '    var.'
    if (alt < 0.) salt = '   undef'
    if ((mszx /= msz) .or. (mszy /= msz)) call abendm('illegal mesh-size')
    is = (ixs-ixl)/mszx;  js = (iys-iyl)/mszy
    if ((is*mszx+ixl /= ixs) .or. (js*mszy+iyl /= iys))   &
      call abendm('ranges conflict')
    if ((mx > mxnb) .or. (my > mynb)) 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 == 0.) then
      read(1,'(a)') buf
      do while (buf(1:1) == '#')
        read(1,'(a)') buf
      enddo
      if (buf(13:16) == '    ') then        !-- v2005
        read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
        if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
        nc1 = ncc1
        if (ncc1 < 200) then
          nc1 = ncc1; ncc1 = nc1 + 800
        else
          ncc1 = ncc1 - 200; nc1 = ncc1
        endif
      else                                  !-- v2018
        read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
        nc1 = ncc1
        if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
      endif
      if ((nc1 >= 1) .and. (nc1 <= 62)) then
        iorg1 = 0;  ispa1 = 0;  ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if (nc1 > 60) korg1 = 0
        if (nc1 == 61) iorg1 = +5400
        if (nc1 == 62) iorg1 = -5400
      endif
      if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg) .or.   &
          (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
        call abendm('Header1 inconsistent')
      read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
      if ((ixs1 /= ixs) .or. (iys1 /= iys) .or. (mszx1 /= mszx) .or.   &
          (mszy1 /= mszy) .or. (mx1 /= mx) .or. (my1 /= my))   &
        call abendm('Header2 inconsistent')
      write(anul,'(f10.3)') vnul1
      if ((anul(9:10) /= '00') .or. (anul(1:1) /= ' '))   &
        call abendm('unable to handle vnul1 value')
      read(1,*) ((t(i,j),i=1,mx),j=1,my)
    else if (alt < 0.) then
      call abendm('Undefined Altitude')
    endif
    close(1)
    write(6,'(a,a)') 'Input filename : ', fnam1(1:lrtrim(fnam1))
    write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
    write(6,'(4x,a)')   &
      '  (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,'(a)') fnam1
    write(9,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
    write(9,*) xs, ys, mszx, mszy, mx, my, alt
    do j=1,my
      jp = js + j + nb
      if ((jp <= 0) .or. (jp > mynb)) cycle
      do i=1,mx
        if (r(i,j) == vnul) cycle
        ip = is + i + nb
        if ((ip <= 0) .or. (ip > mxnb)) cycle
        g(ip,jp) = r(i,j)
        if (alt /= 0.) then
          h(ip,jp) = alt
        else if (t(i,j) /= 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
      enddo
    enddo
    if (nsrc == 1) then
      f(:,:) = g(:,:);  a(:,:) = h(:,:);  g(:,:) = vnuln;  h(:,:) = vnuln
    else
      call dpcini('... Gmerge processing : ')
      do j=1,mynb
        lf = 0;  lg = 0
        do i=1,mxnb
          ndd = nddx
          if (lf < 0) then
            do jm=-nb,+nb
              if ((i+nb) > mxnb) cycle
              if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
              if (f(i+nb,j+jm) /= vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp < ndd) ndd = nddp
              endif
            enddo
            if (ndd /= nddx) then
              lf = 0;  df = -sqrt(float(ndd)) / dnorm
            endif
          else if (lf > 0) then
            do jm=-nb,+nb
              if ((i+nb) > mxnb) cycle
              if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
              if (f(i+nb,j+jm) == vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp < ndd) ndd = nddp
              endif
            enddo
            if (ndd /= nddx) then
              lf = 0;  df = sqrt(float(ndd-1)) / dnorm
            endif
          else
            if (f(i,j) == vnuln) then
              do jm=-nb,+nb;  do im=-nb,+nb
                if (((i+im) <= 0) .or. ((i+im) > mxnb)) cycle
                if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
                if (f(i+im,j+jm) /= vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp < ndd) ndd = nddp
                endif
              enddo; enddo
              if (ndd == nddx) then
                lf = -1;  df = -2.
              else
                df = -sqrt(float(ndd)) / dnorm
              endif
            else
              do jm=-nb,+nb;  do im=-nb,+nb
                if (((i+im) <= 0) .or. ((i+im) > mxnb)) cycle
                if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
                if (f(i+im,j+jm) == vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp < ndd) ndd = nddp
                endif
              enddo; enddo
              if (ndd == nddx) then
                lf = +1;  df = +2.
              else
                df = sqrt(float(ndd-1)) / dnorm
              endif
            endif
          endif
          ndd = nddx
          if (lg < 0) then
            do jm=-nb,+nb
              if ((i+nb) > mxnb) cycle
              if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
              if (g(i+nb,j+jm) /= vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp < ndd) ndd = nddp
              endif
            enddo
            if (ndd /= nddx) then
              lg = 0;  dg = -sqrt(float(ndd)) / dnorm
            endif
          else if (lg > 0) then
            do jm=-nb,+nb
              if ((i+nb) > mxnb) cycle
              if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
              if (g(i+nb,j+jm) == vnuln) then
                nddp = nb*nb + jm*jm
                if (nddp < ndd) ndd = nddp
              endif
            enddo
            if (ndd /= nddx) then
              lg = 0;  dg = sqrt(float(ndd-1)) / dnorm
            endif
          else
            if (g(i,j) == vnuln) then
              do jm=-nb,+nb;  do im=-nb,+nb
                if (((i+im) <= 0) .or. ((i+im) > mxnb)) cycle
                if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
                if (g(i+im,j+jm) /= vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp < ndd) ndd = nddp
                endif
              enddo; enddo
              if (ndd == nddx) then
                lg = -1;  dg = -2.
              else
                dg = -sqrt(float(ndd)) / dnorm
              endif
            else
              do jm=-nb,+nb;  do im=-nb,+nb
                if (((i+im) <= 0) .or. ((i+im) > mxnb)) cycle
                if (((j+jm) <= 0) .or. ((j+jm) > mynb)) cycle
                if (g(i+im,j+jm) == vnuln) then
                  nddp = im*im + jm*jm
                  if (nddp < ndd) ndd = nddp
                endif
              enddo; enddo
              if (ndd == nddx) then
                lg = +1;  dg = +2.
              else
                dg = sqrt(float(ndd-1)) / dnorm
              endif
            endif
          endif
          if (df < 0.) then
            if (dg < 0.) then
              if ((1.+df+dg) > 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 < 1.) then
              if ((1.+df-dg) > 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 < 1.) then
            if (dg < 0.) then
              if ((1.+dg-df) > 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 < 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 < 0.) then
              r(i,j) = f(i,j);  t(i,j) = a(i,j)
            else if (dg < 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)
        enddo
      enddo
      f(:,:) = r(:,:);  a(:,:) = t(:,:);  g(:,:) = vnuln;  h(:,:) = vnuln
    endif
  enddo

  call clspin()
  endfile(9)
  rewind(9)
  write(2,'(a)') '# prog.GMERGE output'
  write(2,'(a8,4x,i4,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
  write(2,'(2i12,4i6,1x,f7.1, 1x,f7.0)') ixl,iyl, msz,msz, mxn,myn, vnuln, 0.
  do j=1,myn
    write(2,'((f7.1,9(1x,f7.1)))') (f(i+nb,j+nb),i=1,mxn)
  enddo
  write(2,'(a8,4x,i4,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn
  write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixl,iyl, msz,msz, mxn,myn, vnuln, -1.
  do j=1,myn
    write(2,'((f7.1,9(1x,f7.1)))') (a(i+nb,j+nb),i=1,mxn)
  enddo
  close(2)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'GMERGE : merge multiple grid data'
    write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km'
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', narea, 'Coord.', nccn, iorgn, korgn, ispan, ispbn
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'gen.new grid', xl, yl, msz, msz, mxn, myn, '    var.'
    write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
    do
      read(9,'(a)',iostat=IER) fnam1
      if (IER < 0) exit
      read(9,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
      read(9,*) xs, ys, mszx, mszy, mx, my, alt
      write(salt,'(1x,f7.0)') alt
      if (alt == 0.) salt = '    var.'
      write(99,'(a,a)') 'Input filename : ', fnam1(1:lrtrim(fnam1))
      write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
        '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
      write(99,'(4x,a)')   &
        '  (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
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end

!-------
  real function avrg(f, i,j, mx,my, vnul, d)
  implicit none
  integer,intent(in) :: i, j, mx, my
  real,intent(in) :: f(mx,my)
  real,intent(in) :: vnul, d
  integer :: n, k, im, jm
  real :: dd, ddp, s

  dd = d * d;  n = nint(d);  s = 0.;  k = 0
  do jm=-n,+n
    if (((j+jm) <= 0) .or. ((j+jm) > my)) cycle
    do im=-n,+n
      if (((i+im) <= 0) .or. ((i+im) > mx)) cycle
      if (f(i+im,j+jm) == vnul) cycle
      ddp = float(im*im + jm*jm)
      if (ddp > dd) cycle
      s = s + f(i+im,j+jm);  k = k + 1
    enddo
  enddo
  if (k == 0) call abendm('Error in function subprogram (avrg)')
  avrg = s / float(k)
  return
  end function