!--------( SELDB : gen. new grid-data from AMDB grid file )--------
  implicit none
  integer :: lwkdir, lrtrim
  character(*),parameter :: AMDBDIR = '/pub/AMDB/DATA/'
  real,allocatable :: f(:,:), g(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1, narea, salt;  character(1) :: ch
  integer :: ldr, l, i, j, n, ip, jp, ix, iy , ksel
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my, ixt, iyt
  integer :: ixl, iyl, mszxn, mszyn, mxn, myn, ixh, iyh
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  real :: vnul, alt, smx, smy, xs, ys, xl, yl, vnul1, dx, dy, f1, f2

  call premsg('SELDB : gen. new grid-data from AMDB grid file')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam2 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call premsg('')
  call gparmi('Select   0 (Survey DB) or 1 (Composite) ==> ', ksel)
  if ((ksel < 0) .or. (ksel > 1)) call abendm('invalid value')
  if (ksel == 0) then
    call premsg(' [Survey DB files :  (areacode).gd ]')
    call system('ls -T0 -CF ' // AMDBDIR // 'grd >&2')
    call gparma('  Enter areacode ==> ', 8, area)
    l = lrtrim(area)
    fnam1 = AMDBDIR // 'grd/' // area(1:l) // '.gd'
  else
    call premsg('[Composite]  a:Hokkaido, b:Kanto, c:Kinki, '   &
              // 'd:Kyushu, e:Ryukyu ')
    call gparma('    Select from  a - e  ==> ', 1, ch)
    if (index('abcde',ch) == 0) call abendm('invalid value')
    fnam1 = AMDBDIR // 'reg/' // ch // '.gd'
  endif
  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
  read(1,'(2i12,4i6,2f8.1)') ixs,iys, mszx,mszy, mx,my, vnul, alt
  write(salt,'(1x,f7.0)') alt
  if (alt == 0.) salt = '    var.'
  if (alt < 0.) salt = '   undef'
  rewind(1)

  allocate(f(mx,my))

  ixt = ixs + (mx-1)*mszx
  iyt = iys + (my-1)*mszy
  smx = float(mszx)
  smy = float(mszy)
  xs = float(ixs)/1000.
  ys = float(iys)/1000.
  write(buf,'(f9.2)') vnul
  if ((buf(9:9) /= '0') .or. (buf(1:1) /= ' '))   &
    call abendm('unable to handle vnul value')
  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6,a8)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'alt '
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6,a8)')   &
    '--- source data ', xs, ys, mszx, mszy, mx, my, salt
  call premsg(out)
  call premsg('')
  call gparma('Enter Name-label for new grid ==> ', 8, narea)
  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 premsg('  Mesh size and numbers to N ?')
  call gparmi('    Unit Mesh size in m :  ', mszxn)
  if (mszxn <= 0) call abendm('invalid value')
  call gparmi('    Number of Mesh to N :  ', mxn)
  ixh = ixl + (mxn-1)*mszxn
  call premsg('  Mesh size and numbers to E ?')
  call gparmi('    Unit Mesh size in m :  ', mszyn)
  if (mszyn <= 0) call abendm('invalid value')
  call gparmi('    Number of Mesh to E :  ', myn)
  iyh = iyl + (myn-1)*mszyn

  allocate(g(mxn,myn))

  call gparma('Enter new grid output filename ==> ', 81-ldr, fnam2(ldr:80))
  call clspin()
  open(2,file=fnam2,status='new')
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'SELDB : generate grid-data for selected area'
  write(6,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', narea, 'Coord.', ncc, iorg,korg, ispa,ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'source grid ', xs, ys, mszx, mszy, mx, my, salt
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'gen.new grid', xl, yl, mszxn, mszyn, mxn, myn, salt

  n = 0
  do while (n < 2)
    read(1,'(a)') buf
    do while (buf(1:1) == '#')
      write(2,'(a)') buf(1:lrtrim(buf))
      read(1,'(a)') buf
    enddo
    n = n + 1
    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)   &
        .or. (vnul1 /= vnul)) call abendm('Header2 inconsistent')
    write(2,'(2a)') '# Source : ', fnam1(1:lrtrim(fnam1))
    write(2,'(a)')  '# prog.SELDB applied to above file'
    write(2,'(a8,4x,i4,4i8)') narea, ncc, iorg,korg, ispa,ispb
    if (n == 2) alt = -1.
    write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
      ixl,iyl, mszxn,mszyn, mxn,myn, vnul, alt
    read(1,*) ((f(i,j),i=1,mx),j=1,my)

    do j=1,myn
      iy = iyl + (j-1)*mszyn
      jp = (iy-iys)/mszy
      dy = float(iy-iys-jp*mszy)/smy
      do i=1,mxn
        ix = ixl + (i-1)*mszxn
        ip = (ix-ixs)/mszx
        dx = float(ix-ixs-ip*mszx)/smx
        if ((ix < ixs) .or. (ix > ixt) .or. (iy < iys) .or. (iy > iyt)) then
          g(i,j) = vnul
        else if (f(ip+1,jp+1) == vnul) then
          g(i,j) = vnul
        else if (dy == 0.) then
          if (dx == 0.) then
            g(i,j) = f(ip+1,jp+1)
          else if (f(ip+2,jp+1) == vnul) then
            g(i,j) = vnul
          else
            g(i,j) = f(ip+1,jp+1) + (f(ip+2,jp+1)-f(ip+1,jp+1))*dx
          endif
        else if (f(ip+1,jp+2) == vnul) then
          g(i,j) = vnul
        else
          if (dx == 0.) then
            g(i,j) = f(ip+1,jp+1) + (f(ip+1,jp+2)-f(ip+1,jp+1))*dy
          else if ((f(ip+2,jp+1) == vnul) .or. (f(ip+2,jp+2) == vnul)) then
            g(i,j) = vnul
          else
            f1 = f(ip+1,jp+1) + (f(ip+2,jp+1)-f(ip+1,jp+1))*dx
            f2 = f(ip+1,jp+2) + (f(ip+2,jp+2)-f(ip+1,jp+2))*dx
            g(i,j) = f1 + (f2-f1)*dy
          endif
        endif
      enddo
    enddo
    do j=1,myn
      write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn)
    enddo
    if (alt /= 0.) n = 2
  enddo
  close(1)
  close(2)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a,$)') 'SELDB :'
    write(99,'(a)')   ' gen. new grid-data from AMDB grid file'
    write(99,'(a,a)') ' Source infile   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', narea, 'Coord.', ncc, iorg,korg, ispa,ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'source grid ', xs, ys, mszx, mszy, mx, my, salt
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'gen.new grid', xl, yl, mszxn, mszyn, mxn, myn, salt
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end