c--------( GALTF : gen. grid-data of flying altitude )--------
      parameter (LXN=1201, LYN=1201)
c
      dimension f(LXN,LYN), w(LXN,LYN), m(LXN,LYN)
      equivalence (w(1,1),m(1,1))
      character*80 wdr, flog, fnam1, fnam2
      character area*8, buf*80
c
      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')
c
      call gparmf('  Specify effecting radius in km ==> ', sr)
      srr = sr * sr
      if(sr.le.0.) call abendm('invalid value')
      call gparma('Enter Name-label for grid-data   ==> ', 8, area)
      call premsg('Specify gridding parameters')
      call gparmi('  Coordinate Number : ', ncc)
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      ispa = 0
      ispb = 0
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      else
        call gparmi2('    Origin Lat. & Long. in minutes : ',
     *               iorg, korg)
        if(nc.eq.72) then
          call gparmi2('    Latitudes of two std.parallels : ',
     *                 ispa, ispb)
		endif
      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.le.0) call abendm('invalid value')
      smx = float(mszx) / 1000.
      ir = nint(sr/smx)
      call gparmi('    Number of Mesh to N :  ', mx)
      if(mx.gt.LXN) call abendm('too large array size')
      if(mx.le.0) call abendm('invalid value')
      call premsg('  Mesh size and numbers to E ?')
      call gparmi('    Unit Mesh size in m :  ', mszy)
      if(mszy.le.0) call abendm('invalid value')
      smy = float(mszy) / 1000.
      jr = nint(sr/smy)
      call gparmi('    Number of Mesh to E :  ', my)
      if(my.gt.LYN) call abendm('too large array size')
      if(my.le.0) call abendm('invalid value')
c
      call gparma('Enter grid-data output filename ==> ',
     *            81-ldr, fnam2(ldr:80))
      call clspin()
      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,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xl, yl, mszx, mszy, mx, my, '    --- '
c
      call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
      do 100 j=1,my
      do 100 i=1,mx
        f(i,j) = 0.
        w(i,j) = 0.
  100 continue
      nls = 0
    1 read(1,'(a)',end=2) buf
      if(buf(1:1).eq.'&'.or.buf(1:1).eq.'%'.or.buf(1:1).eq.'#') then
        call premsg(buf(1:lrtrim(buf)))
        if(buf(1:1).ne.'#') nls = nls + 1
        goto 1
      endif
      read(buf,'(f10.4,1x,f11.4,1x,f8.1)') fi, fk, alt
      call cviken(fi, fk, ye, xn)
      ip = nint((xn-xl)/smx)
      jp = nint((ye-yl)/smy)
      xc = xn-xl - float(ip)*smx
      yc = ye-yl - float(jp)*smy
      do 200 j=-jr,+jr
        if((jp+j).lt.0.or.(jp+j).ge.my) goto 200
        y = float(j)*smy - yc
        do 220 i=-ir,+ir
          if((ip+i).lt.0.or.(ip+i).ge.mx) goto 220
          x = float(i)*smx - xc
          rr = x*x + y*y
          if(rr.eq.0.) then
            wt = 10000.
          else
            wt = srr/rr - 1.
            if(wt.gt.9999.) wt = 9999.
          endif
          if (wt.gt.1000.) then
            if(w(ip+i+1,jp+j+1).ge.(-wt)) then
              w(ip+i+1,jp+j+1) = -wt
              f(ip+i+1,jp+j+1) = alt
            endif
          else if (wt.gt.0.) then
            if(w(ip+i+1,jp+j+1).ge.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
  220   continue
  200 continue
      goto 1
    2 rewind(1)
      do 300 j=1,my
      do 300 i=1,mx
        if(w(i,j).eq.0.) then
          f(i,j) = 9999.9
        else if(w(i,j).gt.0.) then
          f(i,j) = f(i,j) / w(i,j)
        endif
        m(i,j) = 0
  300 continue
c
      call prompt('... gen. flying alt. data : ')
      call dpcent(0, 0)
      nl = 0
    3 read(1,'(a)',end=4) buf
      if(buf(1:1).eq.'#') goto 3
      if(buf(1:1).eq.'&'.or.buf(1:1).eq.'%') then
        call dpcent(nl, nls)
        nl = nl + 1
        goto 3
      endif
      read(buf,'(f10.4,1x,f11.4,1x,f8.1)') fi, fk, alt
      call cviken(fi, fk, ye, xn)
      ip = nint((xn-xl)/smx)
      jp = nint((ye-yl)/smy)
      xc = xn-xl - float(ip)*smx
      yc = ye-yl - float(jp)*smy
      do 400 j=-jr,+jr
        if((jp+j).lt.0.or.(jp+j).ge.my) goto 400
        y = float(j)*smy - yc
        do 440 i=-ir,+ir
          if((ip+i).lt.0.or.(ip+i).ge.mx) goto 440
          x = float(i)*smx - xc
          rr = x*x + y*y
          if(rr.ge.srr) goto 440
          mij = m(ip+i+1,jp+j+1)
          if(i.le.0.and.j.le.0.and.mod(mij,2).eq.0) mij = mij + 1
          if(i.le.0.and.j.ge.0.and.mod(mij/2,2).eq.0) mij = mij + 2
          if(i.ge.0.and.j.le.0.and.mod(mij/4,2).eq.0) mij = mij + 4
          if(i.ge.0.and.j.ge.0.and.mod(mij/8,2).eq.0) mij = mij + 8
          m(ip+i+1,jp+j+1) = mij
  440   continue
  400 continue
      goto 3
    4 close(1)
      call dpcent(nls, nls)
      mij = 0
      do 500 j=1,my
      do 500 i=1,mx
        if(m(i,j).ne.15) f(i,j) = 9999.9
        if(f(i,j).eq.9999.9) goto 500
        if(mij.eq.0) then
          alth = f(i,j)
          altl = f(i,j)
          mij = 1
        else
          if(f(i,j).gt.alth) alth = f(i,j)
          if(f(i,j).lt.altl) altl = f(i,j)
        endif
  500 continue
c
      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,i4,4x,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 600 j=1,my
        write(2,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  600 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,$)') '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,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xl, yl, mszx, mszy, mx, my, '    var.'
        write(99,2000) altl, alth
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'gen.alt grid', 2f10.2, 2i8, 2i6, a8)
 2000 format(' lowest and highest altitudes (m) : ', 2(1x, f7.1))
      end
