!--------( GALTF : gen. grid-data of flying altitude )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), w(:,:)
  integer,allocatable :: m(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, buf
  character(20) :: sdtm;  character(8) :: area
  integer :: ldr, i, j, ir, jr, nls, ip, jp, nl, mij, IER, i1, i2
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ixl, iyl, mszx, mszy, mx, my
  real :: sr, srr, alt, alth, altl
  real :: xl, yl, smx, smy, xn, ye, xc, yc, x, y, rr, wt
  real(8) :: dfi, dfk, dxn, dye

  call premsg('GALTF : gen. grid-data of flying altitude')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input StdLIN data filename ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')

  call gparmf('  Specify effecting radius in km ==> ', sr)
  srr = sr * sr
  if (sr <= 0.) call abendm('invalid value')
  call gparma('Enter Name-label for grid-data   ==> ', 8, area)
  call premsg('Specify gridding parameters')
  call gparmi('  Coordinate Number (+800 for Old-TOKYO) ==> ', ncc)
  nc = ncc
  if (ncc >= 800) nc = ncc - 800
  ispa = 0;  ispb = 0
  if ((nc >= 1) .and. (nc <= 62)) then
    iorg = 0
    korg = (nc-30)*360 - 180
    if (nc > 60) korg = 0
    if (nc == 61) iorg = +5400
    if (nc == 62) iorg = -5400
  else if (nc == 65) then
    iorg = 0
    call gparmi2('  Central Longitude ?  (Deg. Min.) : ', i1,i2)
    korg = i1*60 + i2
  else
    if (nc == 72) call abendm('nc=72 not implemented')
    call gparmi2('  Origin Latitude ?   (Deg. Min.) : ', i1,i2)
    iorg = i1*60 + i2
    call gparmi2('  Origin Longitude ?  (Deg. Min.) : ', i1,i2)
    korg = i1*60 + i2
  endif
  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 :  ', mszx)
  if (mszx <= 0) call abendm('invalid value')
  smx = float(mszx) / 1000.
  ir = nint(sr/smx)
  call gparmi('    Number of Mesh to N :  ', mx)
  if (mx <= 0) call abendm('invalid value')
  call premsg('  Mesh size and numbers to E ?')
  call gparmi('    Unit Mesh size in m :  ', mszy)
  if (mszy <= 0) call abendm('invalid value')
  smy = float(mszy) / 1000.
  jr = nint(sr/smy)
  call gparmi('    Number of Mesh to E :  ', my)
  if (my <= 0) call abendm('invalid value')

  allocate(f(mx,my));  allocate(w(mx,my));  allocate(m(mx,my))

  call gparma('Enter grid-data output filename ==> ', 81-ldr, fnam2(ldr:80))
  call clspin()
  call strdtm(sdtm)
  open(2,file=fnam2,status='new')
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'GALTF : gen. grid-data of flying altitude'
  write(6,'(a,a)') ' Source StdLIN file   : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' Gen. F.Alt grid file : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,f6.2)') '   effecting radius (km) :', sr
  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my '
  write(6,'(4x,a12,2f10.2,2i8,2i6)')   &
    'gen.alt grid', xl, yl, mszx, mszy, mx, my
  call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
  do j=1,my; do i=1,mx
    f(i,j) = 0.;  w(i,j) = 0.
  enddo; enddo
  nls = 0
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%') .or. (buf(1:1) == '#')) then
      if (buf(1:1) /= '#') nls = nls + 1
      call premsg(buf(1:lrtrim(buf)));  cycle
    endif
    i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
    read(buf,*) dfi, dfk, alt
    call cvdiken(dfi, dfk, dye, dxn)
    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
    do j=-jr,+jr
      if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
      y = float(j)*smy - yc
      do i=-ir,+ir
        if (((ip+i) < 0) .or. ((ip+i) >= mx)) cycle
        x = float(i)*smx - xc
        rr = x*x + y*y
        if (rr == 0.) then
          wt = 10000.
        else
          wt = srr/rr - 1.
          if (wt > 9999.) wt = 9999.
        endif
        if (wt > 1000.) then
          if (w(ip+i+1,jp+j+1) >= (-wt)) then
            w(ip+i+1,jp+j+1) = -wt
            f(ip+i+1,jp+j+1) = alt
          endif
        else if (wt > 0.) then
          if (w(ip+i+1,jp+j+1) >= 0.) then
            f(ip+i+1,jp+j+1) = f(ip+i+1,jp+j+1) + alt*wt
            w(ip+i+1,jp+j+1) = w(ip+i+1,jp+j+1) + wt
          endif
        endif
      enddo
    enddo
  enddo
  rewind(1)
  do j=1,my; do i=1,mx
    if (w(i,j) == 0.) then
      f(i,j) = 9999.9
    else if (w(i,j) > 0.) then
      f(i,j) = f(i,j) / w(i,j)
    endif
    m(i,j) = 0
  enddo; enddo

  call dpcini('... gen. flying alt. data : ')
  nl = 0
  do
    read(1,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if (buf(1:1) == '#') cycle
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) then
      call dpcent(nl, nls);  nl = nl + 1;  cycle
    endif
    i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'m');  if (i /= 0) buf(i:i) = ' '
    read(buf,*) dfi, dfk, alt
    call cvdiken(dfi, dfk, dye, dxn)
    xn = sngl(dxn);  ip = nint((xn-xl)/smx);  xc = xn-xl - float(ip)*smx
    ye = sngl(dye);  jp = nint((ye-yl)/smy);  yc = ye-yl - float(jp)*smy
    do j=-jr,+jr
      if (((jp+j) < 0) .or. ((jp+j) >= my)) cycle
      y = float(j)*smy - yc
      do i=-ir,+ir
        if (((ip+i) < 0) .or. ((ip+i) >= mx)) cycle
        x = float(i)*smx - xc
        rr = x*x + y*y
        if (rr >= srr) cycle
        mij = m(ip+i+1,jp+j+1)
        if ((i <= 0) .and. (j <= 0) .and. (mod(mij,2) == 0)) mij = mij + 1
        if ((i <= 0) .and. (j >= 0) .and. (mod(mij/2,2) == 0)) mij = mij + 2
        if ((i >= 0) .and. (j <= 0) .and. (mod(mij/4,2) == 0)) mij = mij + 4
        if ((i >= 0) .and. (j >= 0) .and. (mod(mij/8,2) == 0)) mij = mij + 8
        m(ip+i+1,jp+j+1) = mij
      enddo
    enddo
  enddo
  close(1)
  call dpcent(nls, nls)

  mij = 0
  do j=1,my; do i=1,mx
    if (m(i,j) /= 15) f(i,j) = 9999.9
    if (f(i,j) == 9999.9) cycle
    if (mij == 0) then
      alth = f(i,j);  altl = f(i,j);  mij = 1
    else
      if (f(i,j) > alth) alth = f(i,j)
      if (f(i,j) < altl) altl = f(i,j)
    endif
  enddo; enddo

  write(2,'(a)') '# Flying altitude surface Generated by GALTF'
  write(2,'(a,2(1x,f7.1))')   &
    '# lowest and highest altitudes (m) : ', altl, alth
  write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(2,'(2i12,4i6,1x,f7.1, 1x,f7.0)')   &
    ixl,iyl, mszx,mszy, mx,my, 9999.9, -1.
  do j=1,my
    write(2,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo
  close(2)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a,$)') 'GALTF :'
    write(99,'(a)')   ' gen. grid-data of flying altitude'
    write(99,'(a,a)') ' Source StdLIN file   : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Gen. F.Alt grid file : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,f6.2)') '   effecting radius (km) :', sr
    write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my '
    write(99,'(4x,a12,2f10.2,2i8,2i6)')   &
      'gen.alt grid', xl, yl, mszx, mszy, mx, my
    write(99,'(a, 2(1x,f7.1))')   &
      ' lowest and highest altitudes (m) : ', altl, alth
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end